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ABSTRACT 

The  need  for  tools  capable  of  handling  non-stationariiies  in  the  spectral  content  of 
the  data  has  been  recognized  as  early  as  1946.  The  Wigner-Ville  Distribution  (WD)  has 
been  extensively  used  since  its  introduction  in  1948,  but  suffers  from  some  associated 
problems  (e.g.,  spectral  cross-terms  and  requiring  the  use  of  analytic  signals).  An  alter- 
native Distribution  is  proposed,  which  has  its  origin  in  the  definition  proposed  by  Page 
of  "Instantaneous  Power  Spectrum"  (IPS).  Its  characteristics  are  examined  and,  when 
pertinent,  compared  to  the  WD.  It  is  shown  to  be  less  sensitive  to  the  problems  aflicting 
the  WD,  but  provides  less  frequency  resolution.  The  usefulness  of  a  parametric  (AR) 
version  is  investigated.  Some  typical  test  signals  are  examined,  to  demonstrate  the  per- 
formance and  trade-ofi's  of  IPS  and  its  parametric  version. 
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I.     INTRODUCTION 

A.     THE  PROBLEM 

The  use  of  the  Fourier  Transform  as  a  spectral  description  of  signals  is  a  concept 
whose  usefulness  is  restricted  to  the  class  of  stationar>'  signals.  Though  mathematically 
elegant  and  convenient,  Fourier  decomposition  of  a  signal  can  often  mask  the  true 
spectrum,  since  the  assumed  basis-functions  implicitly  oppose  any  notion  of  time- 
dependency. 

When  the  spectral  content  of  the  signal  changes  with  time,  as  is  often  the  case  in 
fields  such  as  Communications,  Seismology  or  Speech  Processing,  more  powerful  tools 
are  needed.  The  time-dependency  of  the  spectral  content  should  be  apparent  and  meas- 
urable through  the  use  of  a  more  general  type  of  representation  of  the  signal.  Gabor 
(Ref  1]  proposed  such  a  representation  by  introducing  the  concept  of  a  spectral  de- 
scription in  a  joint  time-frequency  plane. 

Since  then,  several  attempts  have  been  made  [Refs.  2.  3,  4,  5]  to  derive  a,  or  the, 
function  capable  of  correctly  describing  the  distribution  of  the  signal's  energy  in  this 
plane.  Short  Time  Fourier  Analysis  (STFA)  has  been  widely  used  and  is  regarded  as  a 
very  convenient  approximation  to  the  "true"  distribution. 

However,  obtaining  this  distribution  is  an  under-determined  problem,  since  an  in- 
finity of  valid  solutions  is  introduced  by  allowing  non-stationary  descriptions.  For 
example,  any  time  signal  can  be  represented  as  a  time-var\'ing  DC  component.  To 
illustrate,  consider  the  following  discrete  time  signal: 

x{n)  =  co<il^  nj  .  (1) 

Two  equally  valid  representations,  in  the  sense  that  each  one  is  capable  of  generat- 
ing the  observed  signal,  are  shown  in  Figures  I  and  2. 

Each  one  of  these  descriptions  explains  the  existence  of  a  sinusoidal  time-series. 
Though  there  is  no  conceptual  distinction  between  them,  the  algorithm  of  analysis  must 
have  the  ability  to  decide  which  one  will  be  taken  as  the  "true"  one.  based  on  an  ex- 
plicitly or  implicitly  built-in  set  of  rules.  The  situation  thus  arises,  where  "one  will  see 
what  one  wants  to  see". 


Figure  1.       x{n)  represented  as  a  constant  sinusoid 
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Figure  2.       .x{n)  represented  as  a  time-varying  DC  component 


Further  complications  appear  when  we  consider  the  fact  that,  given  a  time-series, 
we  cannot  determine  if  the  underlying  process  is  or  is  not  stationar>',  without  using  ar- 
bitrary and  often  inappropriate  assumptions  about  what  is  the  local  behavior.  When  the 
algorithm  of  analysis  produces  a  time-varying  spectrum,  the  fact  remains  that  there  is 
at  least  one  equally  valid  time-invariant  spectral  description  (the  Fourier  Transform 
based),  which  is  arbitrarily  ruled  out. 

B.     APPROACHES 

If  the  signal  is  to  be  represented  in  the  joint  time-frequency  plane  in  a  sensible 
manner,  the  distribution  must  have  some  key  properties.  For  example,  a  shift  in  the 
time-series  should  always  imply  a  corresponding  shift  of  the  spectral  representation 
along  the  time-axis.  Also,  a  multiplication  of  the  time-series  by  a  complex  exponential 


should  result  into  a  shift  of  the  spectral  representation  along  the  frequency  axis.  Without 
these  properties,  physical  interpretation  of  the  representation  can  be  an  impossible  goal. 
This  is  one  of  the  reasons  why  Ambiguity  functions  are  not  used  in  spectral  estimation. 

This  technique  of  imposing  constraints  that  are  felt  needed  in  a  distribution  can  be 
carried  further,  and  the  distribution  (or  one  of  the  distributions)  formed  in  this  refine- 
ment process  defined  as  the  true  distribution.  One  of  the  most  successful  representations 
obtained  by  this  approach  is  the  Wigncr-Ville  Distribution  (WD),  which  has  extensively 
been  used  since  its  introduction  in  1948  [Ref  2]. 

A  different  approach  to  solve  the  indeterminacy  is  to  define,  a  priori  and  unambig- 
uously, what  the  considered  true  distribution  is.  Though  this  approach  lacks  control  of 
the  properties  of  the  resulting  distribution,  it  has  the  advantage  of  more  closely  pre- 
serving physical  meaning  in  whatever  results  it  produces.  A  typical  example  is  the 
definition  of  "Instantaneous  Power  Spectrum"  proposed  by  Page  [Ref  5],  uniquely 
deterniining  a  resulting  distribution  which  is  amenable  to  physical  interpretation. 

C.     OBJECTIVES 

The  Wigner-Ville  Distribution  has  been  extensively  studied  [Refs.  6,  7,  8.  9  ].  and  its 
characteristics  are  fairly  well  understood.  However,  its  use  is  still  hampered  by  some 
problems.  First,  it  has  uneven  performance  for  different  classes  of  signals.  It  performs 
optimally  for  single-component  linear  F.M  (linear  chirp),  but  has  worse  performance  for 
less  regular  spectral  dynamics.  Also,  when  multi-component  signals  are  present,  the  WD 
creates  artifacts  in  the  spectrum,  lying  mid-way  between  true  components.  These 
artifacts  are.  up  to  some  degree,  recognizable  and  treatable  due  to  an  alternating  sign 
[Ref  6].  but  can  severely  mask  results  when  analyzing  more  complex  signals.  The  use 
of  analytic  signals  is  usually  required,  not  only  to  avoid  the  need  for  sampling  at  twice 
the  Nyquist  rate,  but  also  to  prevent  the  appearance  of  undesirable  artifacts  that  would 
otherwise  be  created  by  the  interference  between  positive  and  negative  frequencies 
[Refs.  6,  7,  10]. 

Though  it  has  not  received  proper  attention  in  the  Hterature,  Page's  Instantaneous 
Power  Spectrum  is.  once  properly  understood  and  treated,  a  practical  alternative  to  es- 
timate the  time-var\'ing  spectrum. 

The  study  and  development  of  such  an  alternative  is  the  main  goal  of  this  thesis. 


II.     INSTANTANEOUS  POWER  SPECTRUM  (IPS) 

A.     DEFINITION 

In  an  attempt  to  accommodate  the  notion  of  instantaneous  frequency  content,  Page 
defined  the  Instantaneous  Power  Spectrum  as  the  derivative  of  a  running  energy  spec- 
trum [Ref.  5]: 


p-{t,  J)  =  j^\  S7{J)  p 


(2) 


where 


rt 


S7(J)  = 


s{T)e-^^''^'dx 


(3) 


That  is,  the  Instantaneous  Power  Spectrum  was  defined  as  being,  at  each  fi^equency, 
the  rate  of  change  of  the  energy  collected  by  a  Fourier  transform  taken  fi"om  —  oo  .  up 
to  the  time  of  analysis.  This  concept  was  later  extended  by  Levin  [Ref.  3],  who  intro- 
duced a  complementar\'  backward  run,  similarly  defined  as: 


-  C    I   f,+,^  ,2 


p''{i^f)  =  -^\st(J) 


(4) 


where 


/*oo 


stiJ) 


5(t)£-^''"^Vt  , 


(5) 


and  defining  the  Instantaneous  Power  Spectrum  as  the  average  of  these  two  runs: 


IPS{t,J)  =  \[p-{t,J)  +  p^{t,f)\. 


(6) 


This  expression  can  be  put  in  the  form  [Refs.  3,  1 1  ] 

IPSU,  J)  =  Real  [5*(0-S(y)e^'"-^'] 


(7) 


where  S(J)  is  the  Fourier  Transform  of  5(/). 

IPS  can  thus  be  seen  to  be  the  real  part  of  Rihaczek's  distribution  [Ref  4],  which 
allows  one  to  think  of  IPS{t,f)  as  being  the  energy  density  entering,  at  time  /,  an  infi- 
nitely narrow  filter  centered  at  frequency/  An  enlightening  treatment  of  this  point  can 
be  found  in  Ackroyd  [Ref  12]. 

1.     Time  Domain 

Let  us  now  consider  the  following;  we  want  a  function  of  /,  t  whose  Fourier 
Transform  is  IPS{i,f)  .  That  is,  we  want  G(/,  t)  such  that 


IPS{t,f)  = 


G{[,  t)(?"^"^''-^Vt 


(8) 


Then,  using  (7) 


G(^T)  =  ^ 


[s{i)S*{f)e-''''^'  4-  s{!fSiJ)e^''^':\  e^^^'df 


—oc 


(9) 


from  where,  with  the  appropriate  variable  substitutions  we  obtain 


oo  oo 


G{i.r)  =  ^ 


sii 


-n-fn. 


s  {[  -  v)<?"'  "^  V;/  +  5  (/) 


*J  ^  k) 

— oo  — oo 


s{i  +  }])e~^'^''^y)] 


e''^''^'df.      (10) 


Interchanging  the  order  of  integration,  we  get 


G(M)  = -^  W/)5  (/ -  t)  +  5  (05(/ +  t)]  . 


(II) 


IPS  can  now  be  restated  as: 


IPS{t,f)  =  -^ 


ls{t)s\t  -  t)  +  s\t)s[t  +  t)]  e-^^^'^'dj 


(12) 


— oo 


2.     Frequency  domain 

Alternatively,  a  dual  expression  can  be  derived,  expressing  IPS  in  terms  of  the 
Fourier  Transform  of  the  signal. 


IPS(l,f)  =  4r 


5(0/(/-t)e"^'^^-^VT  +  Y 


s\i)s{i  +  T)e~^^''^'dx 
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With  the  appropriate  variable  substitutions,  we  get 


IPS{i,f)  =  -^ 


Sif+y)S{f}e^'''''dy^^ 


S{f-y)S(J)e^'''"dy  (14) 


—OG 


and,  finally 


IPSit,f)  =  -^ 


lSif)S*if-  y)  +  S*{/)Sif+  >•)]  e^'^^'^dy 


(15) 


3.     IPS  for  discrete  signals 

The  discrete  version  of  IPS  follows  directly  from  (  12  )  as 


(16) 


IPS[n.  6)  =  ^     >    [s{n)s\n  -  k)  +  s\n)s{n  +  ky]e~^^^. 

A'= — >o 

A  formal  derivation  of  (  16  )  can  be  found  in  Appendix  A. 

B.     RELATIONS  WITH  OTHER  TIME-FREQUENCY  DISTRIBUTIONS 

1.     Rihaczek's  Distribution 

The  Rihaczek  Distribution  was  proposed  in  1968  as  the  true  energy^  distribution 
in  an  attempt  to  unify  existing  results. 

Formally,  the  distribution  is  [Ref  4] 


£(/.  J)  =  s[i)S*{/)e-'-^' 


(17) 


where  s{i)  is  the  analytic  signal  being  analyzed,  and  S{f)  is  its  Fourier  Transform. 

As  was  derived.  £(/„.  f^)  is  the  complex  energy  density  at  the  point  [Iq./o)  in  the 
time-frequency  plane.  A  similar  expression  had  been  considered  by  Levin,  who  defined 
the  Complex  Insianianeous  Power  Spccirum  as  the  complex  conjugate  of  (   17  ). 


Rihaczek  used  the  concept  of  "complex  energy",  whose  real  part  is  the  real  energ\'.  while  the 
imaginap.'  part  is  the  reacti\e  energy. 


Using  Levin's  notation,  the  Complex  Instantaneous  Power  Spectrum  is 


(18) 


where  s{i)  is  the  complex  envelope  of  the  real  signal  under  analysis. 

We  thus  see  from  (  7  )  that,  though  addressing  a  broader  class  of  signals,  IPS 
is  both  the  real  part  of  Rihaczek's  Distribution  and  the  real  part  of  Levin's  Complex 
Instantaneous  Power  Spectrum.   Formally, 


lPS[t,J)  =  Real[£(/,  /)] 


(19) 


2.     Cohen's  Generalized  Phase-Space  Distribution  Functions 

In  1966,  Cohen  introduced  a  generalized  class  of  time-frequency  signal  repres- 
entations, given  by 


oo        oo       oo 

r    r    r 


K^  J)  = 


0(i>,  t)  s{t^  +  -^  )5*(/,  -  y  )  e^^""^"''  -"-'^civdi^dr  (20) 
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where  the  choice  of  cp(u,  t)  will  determine  the  resulting  distribution  [Ref  13]. 
If  we  slightly  rearrange  (  20  )  to  read 


oo        oo 

t*      f* 


c{t,  f)  = 


<^{v,  t) 


%I      %) 

— oo     — oc 


oo 


s{t,+^)s{t,-^)^~''''^dt 


'—V 

— oo 


e-^^''^"^'^dvdx ,       (21) 


we  can  interpret  it  as  the  double  Fourier  Transform   of  T(i),  t)  [Ref  14],  where 


oo 


^{0,  t)  =  (D(u,  t) 


sih  +  y  )-^  (^1 


■^)e^'''^^Vr, 


— oo 


(22) 


We  can  now  recognize  (  22  )  as  a  general  expression  for  the  well  known  Ambi- 
guity function,  which  can  serve  as  a  general  definition  of  the  time; frequency 
autocorrelation  function  [Refs.  3,  4  ]. 


Each  choice  of  weighting  function  0(0,  x)  in  (  20  )  provides  a  different  distrib- 
ution, because  each  one  defines  the  combined  autocorrelation  (  22  j  in  a  somewhat  dif- 
ferent way. 

Let  us  see  which  distributions  result  from  three  different  realizations  of  O(u.t) 
(complex  exponential,  constant,  and  real  sinusoid)  [Refs.  4,  13,  14  ]. 

a.     Complex  exponential 


jlzv- 


0(d,  t)  =  e      2 


(23) 


The  combined  autocorrelation  becomes 


oo 


T(u,  t)  =  ^■"•^^ 


5(/,   +^)5*(/,  -^)^''"''V/, 


5(r, +  — )5  (/i-  — y     ^'     2  V/, 


(24j 


and,  defining  a  new  variable  /  =  /,  +  —.  we  obtain 


Y(l).  t)  = 


5(/,)5  (/,  -T)e'^''''''di 


J 

— oo 


(25) 


The  resulting  distribution  is 


CX5         oo        c» 

r    r    r 


c{r,  J)  = 


s{i{)s\t,  -  t)^'''^'-''  e-^^''^'''^'^di,dxdv  , 


J      J 

— oo     — oo    — oo 


(26) 


c{i,f) 


s[i{)s\t^  -j)e  ^^"^'^6(1^  —  i)dt^dx 


— oo    — oo 


(27) 


— oo 


s{t)s\t  -  T)e~^^'''^dx 

which  is  the  well-known  Rihaczek  Distribution. 

Though  this  definition  of  the  combined  autocorrelation  is  probably  the 
most  intuitive  one.  its  lack  of  mathematical  symmetry  has  the  undesirable  effect  of 
making  the  distribution  complex.  It  has  extensively  been  used  in  radar  theor\\  since 
(  25  )  is  exactly  Woodward's  definition  of  the  Ambiguity  function  [Ref.  3]. 

b.     Constant 


(D(u,  t)  =  1 


(28) 


The  combined  autocorrelation  becomes 


"Viv,  t)  = 


5(r, +-^)5*(/,-^)e^''"^Vf, 


— oo 


(29) 


The  resulting  distribution  is 


oo        oo       oo 


c{i,  J)  = 


5(/,  +  —  )S    (/, 


^^^^•2..;,^-y2.(.r+./)^^^^^^^^  (30) 


— oo     — oo    — oo 
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c(^y) 


5(/,  +  -^  )5*(/,  -  -^  )^-^''^^^d(/,  -  /y/,^T 


— oo   — oo 


(31) 


5(r  +  y)5  (r-  — )£  -^     Vt 


which  is  the  Wigner-Ville  Distribution.  This  definition  of  the  combined  autocorrelation, 
though  only  slightly  different  from  the  one  in  (  25  ),  has  the  needed  symmetr\'  to  make 
this  distribution  real,  and  mathematically  ver\'  convenient  (when  dealing  with  analytic 
signals). 

f.      Cosine  function 


0(u,  t)  =  cos(;:t;T) 


(32) 


The  combined  autocorrelation  becomes 


r 


^'{v,  t)  =  cos(;rur) 


s{i,+^)s\i,-^)e^''''\Jt, 


— oo 


(^--  +  e-^"""-') 


oo 


5(/,+^)5*(/,--^)£^''"^V/, 


— oo 


.(/,+f)/(/,--f)^^-^'-H^ 


(33; 


+ 


s{t,+j-)s\t,-f)e^-"^'^-T'dt, 


11 


we  get 


Substituting  t  for  /,  +  —  in  the  first  integral,  and  /  for  /,  —  —  in  the  second, 


H^(f^,T)  =  -^ 


s{t)s  {t  -  r)e'~''''cit  + 


s\t)s{t  +  T)e^^''''dt 


— oo 


— oo 


(34) 


[5{t)s\t  -  t)  +  s\t)s{t  +  r)y^''"dt . 


The  resultine  distribution  is 


oo        oo       oo 

r    r    r 


cU.  D  = 


-^  [5(/i)5>,  -  t)  +  /(/,Kr,  +  T)y''^% 


2r:v:,-j2-{vl+'J) 


dt^dvdr 


— oo     — oo    — oo 


oo       oc 


\  {s{i,)s\i,  -  t)  +  5*(/,K/,  +  x)Y^^''--^b{t,  -  t)dt,dT  (35) 


-oo    — yo 


[s{i)s\i  -  t)  +  s\t)s{i  +  r)\e'^^'''-^dx 


J 

— oo 


which  is  the  expression  we  had  found  earlier  for  IPS.  In  this  case,  the  definition  of  the 
combined  autocorrelation  also  possesses  the  symmetry  needed  for  a  real  distribution. 
The  kernel  that  generates  IPS  ( cl)(u,  t)  =  cos(;rL)T) )  has  been  considered  by  Cohen 
[Rcf  13].  and  shown  to  generate  the  Margeneau-Hill  distribution  [Ref  15],  well  known 
in  quantum  mechanics  theor\'. 
3.     Wigner-Ville  Distribution 

Since  the  Wigner-Ville  Distribution  (WD)  has  enjoyed  wide  acceptance,  we  will 
define  more  closely  the  relations  between  IPS  and  WD.  As  an  indirect  connection,  we 
can  relate  both  to  the  Rihaczek  Distribution.  From  (  29  )  the  WD  implicitly  defines  the 
combined  autocorrelation  as  [Ref  4] 
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H'(o,t)  = 


^U, +T)-^*(^i-T)^''"''^^i 


— oo 


— /rruT 


5(M/a,  -  ty'^^^V/, 


(36) 


Hence,  the  Wigner-Ville  Distribution  is  given  by  [Ref.  4] 


IVDU,  ./)  =  f. 


-y-ii 


5(/,)/(/, -T^^^Vr 


(37) 


where  F,_[  .  ]  denotes  the  2-D  Fourier  transform.  From  the  convolution  property  of  this 
operator,  and  using  (  25  ), 


r  r 


WD{i,J)  =  F^^le-^'''']**F,^, 


5(/,)5   (/i  -  t)^'^*- V/ 


— oo 


^..^['-'"'"T**^^^) 


_  /■'^"/^  ** ., 


^(^y) 


(3S) 


where    **    stands  for  2-D  convolution,  and  e{[,  f)  is  the  Rihaczek  Distribution.  From 
(  38  ),  using  (  19  )  and  the  realness  property  of  WD  (see  Section  II-C),  we  obtain 


rrz)(/,  y)  =  Real[^■''^■'"'**£(^y)] 

IPS{i,  f)  =  RealL  d{t,  f)  **  c{t,  y)]  . 


(39) 


We  thus  see  how  the  two  distributions  can  be  obtained  from  the  Rihaczek  Distribution 
by  an  appropriate  choice  of  the  function  convolved  with  c{t,  J) . 
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We  found  earlier  that  each  choice  of  T(u,  t)  in  Cohen's  class  of  Distributions 
(  20  )  yields  a  particular  definition  of  the  combined  autocorrelation  function 


^{v,  t)  =  (D(u,  t) 


5(/, +y)5>i 


■f  )e^'''"^dt, 


(40) 


— oo 


and  that  the  2-D  Fourier  transform  of  {  40  )  resulted  in  a  particular  distribution,  de- 
pending on  0{d,  t). 

We  see  now  that  (  40  )  is  the  product  of  two  functions,  <l>{o,  t)  and  the  integral. 
Its  double  Fourier  Transform  is  thus  a  double  convolution  in  the  transform 
domain. [Ref   14] 

That  is: 

c(/,y)  =  F,,,[4^(o,T)] 


=  f... 


cD(u,  t) 


5(r,  +  ^M/,-^y''^^V/i 


— oo 


=  F,,,[cI)(o,r)]**F,,, 


5(/,+-f  )5Vi-t)^''"''^^1 


=  F,,,[cl)(u,T)]**n-Z)(r.y). 

In  the  particular  case  of  IPS, 

IPS{t,  J)  =  F^^  X  cosinor)-]  **  ]VD{t,  J)  . 


Therefore, 


IPS{t,J)  =  cos{Anfi)**WD[t,J), 


(41) 


(42) 


(43) 


which  we  could  have  deduced  directly  from  (  39  )  and  the  fact  that  the  WD  is  always 
real,  as  follows: 


IPS{t,  j)  =  Real[£(r,  J)-] 


=  i"  w^  y)  +  £  (/,  y)] . 


(44) 
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Using  (  38  ), 


IPS{i,  /)  =  -j  [e'^'^""^'  **  lVD{i,  J)  +  e'^''^'  **  IVD{[,  /)] 
=  cos{4nft)  **  \VD{t,  J)  . 


(45) 


At  this  point,  we  would  like  to  be  able  to  express  the  WD  in  terms  of  IPS,  with 
a  relation  of  the  form: 


wD{t,  J)  =  oupsii,  m , 


(46) 


where  0[  .  ]  is  an  arbitrary"  operator.    Unfortunately,  the  inverse  filter  does  not  exist, 
since 


f., 


1 


C0S(7rL)T) 


(47) 


is  not  defined  over  the  required  range  of  u  and  r  . 
4.     Short-Time  Fourier  Transform 

Though  the  Short-Time  Fourier  Transform  (STFT)  can  not  be  considered  a  true 
Time-Frequency  Distribution  [Ref  14].  its  widespread  use  fully  justifies  the  study  of  how 
h  relates  to  IPS.  This  hs  been  addressed  in  [Ref  11],  and  will  be  done  here,  after  which 
some  indirect  results  will  emeree.  Let  us  thus  consider  the  STFT  as  defined  in  (  48  ): 


STFT{t,co)=  \F{{,(o)\ 


s{i  +  t)vv*(t)^''"Vt 


— oo 


oo        oo 


(48) 


s{i  +  r)s*{[  +  d)w  {T)w{u)e  -^^'edrdo 


— oo    — oo 
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Substituting  s{i  +  r)  and  vv*(t)  by  their  Fourier  definitions,  we  have  that 


STFT{i,co)  = 


s  {i  +  v)w{d) 


SiOe^'^'^M 


— oo     — oo 
oo        oo 

r    ^ 


J 

— oo 


lv\p)e-^P'dpe-^'"'e^'^''drdo 


— oo 


s\t  +  d)w{d)S{p  +  a>)^r*(p)/'"+^^V'"Vi;^p 


— oo    — oo 
oo       oo 


s\i  +  v)S{p  +  oj)e'^^"'^^'^'^w{D)W\p)e-^''dvdp 


— oo    — oo 


(49) 


We  thus  see  that  the  STFT  is  a  time-frequency  smoothed  version  of  the  Rihaczek  dis- 
tribution of  the  signal.  The  smoothing  function  is  the  Rihaczek  distribution  of  the 
window  used  to  compute  the  STFT.  Since  the  right-hand  side  of  (  49  )  is  ahvays  real  and 
positive,  and  the  previous  relations  are  valid  for  arbitrary  signals  and  windows,  we  con- 
clude that  windowing  a  Rihaczek  distribution  with  a  Rihaczek  distribution  guarantees  a 
real  and  everwhere  positive  distribution.  This  result  will  be  important  when  addressing 
the  issue  of  positivity. 

C.     BASIC  PROPERTIES  OF  IPS 

As  was  already  mentioned,  and  is  easily  seen  from  (  20  ),  it  is  the  choice  of  a  par- 
ticular Ofi),  t)  that  determines  the  particular  distribution  and,  hence,  its  properties.  It  is 
thus  desirable  to  establish  direct  connections  between  the  properties  of  0(0,  t)  and  the 
properties  of  the  resulting  distribution.  This  issue  has  been  addressed  [Ref  14].  and  we 
will  in  the  most  part  only  state  results,  as  apphcable  to  IPS.  The  basic  properties  of  IPS 
are: 

1 .    Time  shift 


If  co(/-)  =  s[t  —  Q,  then 

2.    Modulation 

If  a)(/)  =  5(/)^-'2''V,  then 


IPSJi,f)  =  IPS,ii-i„J) 


(50) 


IPSJt,J)  =  IPS,{i,f-fo). 


(51) 
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3.   Marginal  in  Time 


IPS{i,  f)  df=  I  sii) 


4.   Marginal  in  Frequency 


oo 


IPS{t,J)di  =  \S{f)\\ 


— oo 


where  S{J)  =  Fls{i)']  . 

5.  Realness 

6.  Instantaneous  Frequency 


JPS{i,  j)  is  real  for  all  /,  /. 


f-IPS{i,  f)df 


%j 

— oo 


S{l) 


=fi{t)      {instantaneous  frequency) 


7.   Group  Delay 


oo 


fIPS{[,J)di 


SiJ) 


=  To(J)      {group  delay^ 


8.   Zero  Power 


s{t^)  =  0  =>  //>5(ro,  y)  =  0 
5(/o)  =  0  =^  lPS{t,  /o)  =  0 


(52) 


(53) 


(54) 


(55) 


(56) 


(57) 


The  Wigner-Ville  Distribution  also  respects  properties  1-7  [Ref.  6],  but  only  a 
weaker  version  of  property  8. 

D.     FURTHER  PROPERTIES 

Though  the  properties  presented  in  Section  II-C  constitute  the  basic  set  of  rules  to 
which  IPS  is  bound,  some  further  analysis  must  be  made.  The  efiects  of  linear  operators 
should  be  investigated,  to  fully  understand  the  behavior  of  IPS. 


~  This  property  is  valid  for  anal\lic  signals  only  (  see  Appendix  B  ). 
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I.     Windo>\ing  in  the  time-domain 

The  issue  here  is  to  determine  how  IPS  is  afiected  by  windowing  the  original 
time  signal. 

We  will  at  this  point  introduce  the  notation: 


ipsj^t,  /)  =  ^  ips^{t,  y)  +  ^  ips^it,  J) 


(58) 


where 


IPSait,  f)  = 


d{t)d\t  -  T)e-^^''^'dT 


— oo 


(59) 


and 


lFSZ{t,  f)  = 


d\t)d{t  +  i:)e~^^''^'d 


(60) 


If  the  sisnal  is  windowed,  that  is: 


d{t)  =  5(/)vv(0  , 


(61) 


then  the  IPS  of  the  windowed  sienal  will  be 


IPSltJ)  =  ^ 


— oo 


— oo 


{_d{i)d\t  -  t)  +  d\i)d[i  -  x)-]e-^^''^'dx 

\_s{t)w{t)s\t  -  T)w\t  -  t)  + 
+  s\t)w\t)s{i  +  t)vv(/  +  t)](?"-^'^"-^Vt  . 


(62) 
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Separating  the  two  terms, 


IPS,{[,  f)  =  -^\      ls{i)s\i  -  T)Xw{i)w\i  -  t)]^-^'''^^Vt 


+  - 


ls\r)s{i  +  t)][w*(Ovv(/  +  t)]^"-''^''-^Vt 


— CX3 


r 


— oo 


F  \iPs;:\F  \iPs::]e-^^''-^'dT 


(63) 


+ 


-1 


F    [IPS^^F    llPS^'\e-^'^''-^'dx 


oo 


/?s;(MO/Ps;(/./-uyu  +  ^ 


/PS;(/,  u)IPSt{iJ-  o)do  , 


— oo 


where  F  [.]  is  the  inverse  Fourier  Transform  operator.   Hence 


ips^i,  yi  =  ^  [iPs;{t.  J)  *  iPSZ{t.  /)]  +  y  [^^'^^^■'(^  7)  *  JP^tU,  f)\  ,     (64) 


vvhere  *  denotes  convolution  in  the  f  variable. 

That  is.  while  for  the  Wigner-Ville.  windowing  the  time  signal  implies  the  con- 
volution along  the  frequency  axis  of  the  signal's  WD  with  the  window's  WD  [Ref  6]. 
for  IPS  we  have  a  sum  of  two  convolutions:  the  first  term  of  the  signal's  IPS  [IPS;)  is 
convolved  with  the  first  term  of  the  window's  IPS  {IPS'),  and  a  similar  convolution  is 
applied  to  (IPS;)  and  {IPS;)  . 

A  special  case  occurs  when  the  window,  the  signal,  or  both,  possess  conjugate 
symmetn.'  about  /,  the  time  at  which  we  are  computing  IPS  (  i.e., 
s{t  +  t)  =  s'{i  —  t),  \v{t  +  t)  =  w\i  —  t),  or  both  ).   If  w{t)  possesses  that  symmetry,  then 


iPSZit,  J)  =  ipstu,  D  =  /PSJ^  J) 


(65 
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and.  by  linearity  of  the  convolution  operator,  (  64  )  becomes 


IPS^t,  J)  = 


IPS,ii,  u)IPS^{i,  f-  o)do 


(66) 


We  thus  see  that  in  this  case  IPS  and  WD  respond  similarly  to  the  windowing 
of  the  data:  the  IPS  of  the  windowed  signal  is  the  convolution  (in  frequency)  of  the  IPS 
of  the  signal  with  the  IPS  of  the  window. 

In  summar}',  the  IPS  of  a  windowed  signal  is  a  frequency-smoothed  version  of 
the  IPS  of  the  unwindowed  signal.  The  time-resolution  is  not  affected. 
2.     Convolution  in  the  Time-domain 

The  question  now  is:  how  is  the  IPS  of  a  signal  affected  if  the  signal  is  pre- 
processed  with  a  filtering  operation?   So,  let 


d{i)  =  s{i)  *  h{t) . 


Then. 


IPS,{!,  f)  =  Real  [[/(/)  *  h\t)-]rW)  *  h[t)-]  e'^-"^'] 
=  Real  [ls\t)  *  h\t)-]S{j)H{J)e^^''^'] 
=  Real  [[5*(/)S(/)£^''"-^']  *  lh\i)H{f)c^^^^'-\\  , 


(67) 


(68) 


where  we  used  the  fact  that 


e^^''\a{t)  *  bit)-]  =  e^^"--^' 


a{r)b{t  —  r)dT 


— oo 


(69) 
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We  thus  see  that  (contrary  to  what  happens  in  the  WD-  )  filtering  the  signal 
does  not.  in  general,  correspond  to  the  convolution  in  time  of  the  signal's  IPS  with  the 
IPS  of  the  impulse  response  of  the  filter.  What  is  observable  from  (  68  )  is  that  the  II\S 
of  the  filtered  signal  is  now  the  real  part  of  the  convolution,  along  the  time  axis,  of  the 
two  Rihaczek  Distributions.  The  net  efTect  produced  on  IPS  is,  thus,  still  a  smoothing 
operation  along  the  time  axis.  The  frequency  resolution  is  not  afTected. 
3.     IPS  and  Moyal's  Formula 

As  is  often  stated  in  the  literature,  IPS  does  not,  in  general,  respect  Moyal's 
formula  [Ref.   16].   When  applied  to  the  WD,  Moyal's  formula  states  [Ref   17]: 


IVD.ii,  J)\VDXt,  J)didf=  <  s,g  >  <  s,s  > 


(70) 


— oo  — oo 


where  <  s,g  >  is  the  inner  product  of  5  and  ff  .   Let  us  now  see  what  happens  in  the  IPS 
case: 


oc       oo 

r     r 


IPSJi,  f)IPS,.(i,  f)didf= 


i-jC;  or.' 

r    ^ 


— oo     — oo 


is{l)s*{i  -  r)  +  /(/),s(/  +  T)']e~^^''^'' dT 


(71) 


J 

— oo 


[jg{i)g\i-v)-^g\t)g{t+v)-\e  ^^'^^Vu 


— oo 


didf. 


'  1  he  WD  of  a  filtered  signal  is  the  convolution  along  the  time-axis  of  the  original  signal's 
WD  with  the  WD  of  the  Impulse  response  of  the  filter  [Ref  6j. 
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oo       oo 


IPS,ir,  niPSJi,  f)dtdf= 


— oo    —CO 


oo        oo        oo 

r    r    r 


[s{i)s  (/  -  t)  +  5  {i)s{t  +  rmg{l)g  U-o)+g  {t)gil  +  u)] 


J      J 

— oo     — oo    — oo 


-J27:J\r+v) 


df 


— oo 


dxdtdv 


(72) 


oo       oo 


[5(/)5  (/  -  t)  +  5  {t)s{t  +  i:)-][_g{t)g  (r  +  t)  +  ^  {t)g{t  -  t)]  dt  dx 


%1        *J 

— oo    — oo 


1  *  1 

—    <s,g>     <  s.g>  -\-—    <  S,g>  <  S.g> 

4  4 


+  y  Real 


s{t)g{t) 


s  {t-  x)g  {t  +  x)dx 


di 


and,  finally 


oo       oo 

r    r 


IPS,(i,  f)IPS.{i,  f)didf-. 


— oo    — oo 


1  *       1 

—  <s,g>   <s,g>    +yReal 


sir)g{t) 


s  {t-  x)g  (/  +  x)dx 


— oo 


dt 


(73) 


The  result  confirms  that  M oval's  formula  is  not,  in  general,  valid  for  IPS.  Only 
in  the  trivial  case  of  one  of  the  functions  being  a  real  constant  would  the  second  term 
in  the  right  hand-side  of  (  73  )  reduce  to  —  <s,g>  <s,g>'',  and  hence  would  IPS  obey 
M oval's  formula. 
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4.     Recovery  of  time-signals 

To  address  the  problem  of  recovering  a  signal  from  its  IPS,  we  will  find  it  con- 
venient to  consider  two  separate  cases. 

a.  Infinite  duration  signals 

For  real  signals  with  5(0)  #  0,  recovery  of  the  signal  is  always  possible  ex- 
cept for  a  sign  indeterminacy,  by  considering 

IPS{i,  0)  =  Real[5(0S(0)] 

(74) 
=  S{0)s{t) 

and,  since  |  S(0)  |  can  be  found  by  using  Property  4  (  53  ),  only  the  sign  of  5(0)  remains 
undetermined.  No  similar  results  are  available  for  complex  signals  of  infinite  duration. 
However,  if  s{t)  possesses  conjugate  symmetry  around  some  /g,  then  we  can  also  recover 
complex  signals  of  infinite  duration  up  to  a  phase  term,  by  considering  P'  [  IPS{i(^.J)  ]. 
Similarly,  if  S{J)  possesses  conjugate  symmetry  around  some  j^  and  5(0)^^0,  we  can 
again  recover  complex  signals  of  infinite  duration  up  to  a  phase  term,  by  considering 
FlIPS{!j,)l 

b.  Finite  duration  signals 

For  finite  duration  signals,  recovery  of  real  or  complex  signals  is  always 
possible  except  for  a  constant  phase  term  (or  a  sign,  in  the  case  of  real  signals). 
Let  us  denote  the  starting  time  by  i  =  0.  From  (  12  ) 

F~\lPS{i),  /j]  =  y  [  s{0)s\  -  t)  +  /(0)5-(t)  ]  (75) 

and,  since  5(  —  t)  =  0  for  t  >  0  , 

F~'  UPS{0,  m  l.>o  =  y  ^*(0)5(t)  (76) 

which,  by  property  3  (  52  ),  lea\es  an  unknown  phase  term. 

E.      IPS  -  A  GENERALIZATION  OF  THE  VMENER-KHINCHIN  THEOREM 

By  redefining  auto-correlation,  IPS  offers  a  very  straightforward  extension  of  the 
Wiener- Khinchin  theorem  to  the  non-stationary  case,  as  we  shall  see. 
Let  us  define  the  auto-correlation  function  of  a  signal  as 

R[i,  r)=\  E{s{t)s\[  -  t)  +  s\[)s{i  +  t)}  .  (77) 
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We  have,  then,  using  (  12  )  and  by  linearity  of  the  convolution  operator  the  follow- 
ing pair  of  relations: 


•     E{IPSii,f)}  = 


R{t,r)  e-J^''^dT 


R{i,  t)  = 


E{IPS{1,  J)]  e^'^f^df 


These  relations  can  be  read  as  a  generalization  of  the  Wiener-Khinchin  theorem  to 
non-stationary  signals,  since  for  stationary'  signals  they  promptly  reduce  to  that  theorem. 
U  s{t)  is  a  stationary"  signal,  then 


R{t,  t)  =  Y  E[s[t)s\i  -  t)  +  s\t)s{i  +  t)} 

=  -^  E{s{t)s\t  -  t)}  +  -|-  E[s\t)s{i  +  t)} 

=  ^R[t)  +  ^R\-t) 
=  R[t) 


(78) 


E{lPS[i.f)]  = 


R[T)e~^^^^'dT  =  I  SiJ) 


— oo 


(79) 


The  expected  value  of  IPS  is,  thus,  for  stationar}'  signals,  the  usual  Power  Spectral 
Density  (  PSD  ). 

Since  in  the  non-stationan.'  case  we  cannot  tr\'  to  infer  ensemble  averages  from  time 
averages,  the  IPS  of  a  signal  is  one  realization  of  the  generalized  Wiener-Khinchin  re- 
lations of  the  underlying  stochastic  process. 

A  different  view  of  IPS  when  applied  to  stochastic  processes  can  be  found  in  Grace 
[Ref  11].  where  E[lPS[i,  J)]  is  shown  to  be  the  Fourier  Transform  of  Loeve's  General- 
ized Power  Spectral  Density. 
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111.     RETHINKING  IPS 

A.     ON  THE  USE  OF  WINDOWS 

IPS,  as  given  in  (  12  ),  sufTers  from  intense  ringing  at  all  frequencies  "touched"  by 
the  signal.  This  same  eOect  is  one  of  the  reasons  why  Rihaczek's  Distribution  did  not 
fmd  much  acceptance.  An  example  can  be  seen  in  Figure  3,  where  a  linear  chirp  was 
used  as  the  test  sisnal. 


Figure  3.       IPS  for  a  linear  chirp 


This  ringing  elTect  can  easily  be  understood  by  slightly  rearranging  (  12  )  to 


IPS{i,J)  =  ^ 


s{i) 


"—      %/ 


s\[  -  t)c~^^'''-^' dx  +  s\i) 


s{i  +  T)e~^^''^'dx 


(80) 
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With  the  appropriate  changes  of  variable,  we  get 


//'S(/,y)  =  -^ 


5(/) 


/(r,)^-^'^-^"''V/,  +  /(0 


— oo 


<hW 


ln/{t-  /j) 


du 


(81) 


=  Real  [  s\t)  [5(0  *  e^^'^-^']  ]  . 

But  s{t)  *  e^^"^'  is  the  output  of  a  filter  whose  impulse  response  is  a  complex  exponential. 
IPS  is,  thus,  for  each  frequency,  the  real  part  of  the  product  of  the  signal's  complex 
conjugate  and  the  output  of  an  infinitely  narrow  non-causal  filter  (non-causal  oscillator) 
(Figure  4).   (Ref   12] 


s(t) 

jarft 

h(t>e 

. — ' 

' 

(x^ 

Real   (  ) 

IPS(tf)^ 

• 

I 

» 

(  ) 

- 

Figure  4.       Model  of  IPS 

The  ringing  is  hence  unavoidable,  and  will  persist  until  the  filter's  impulse  response  be- 
comes negligible.  We  should  consider  stable  filters,  with  decaying  impulse  responses,  if 
we  want  to  diminish  this  efiect.  This  reasoning  leaves  us  very  close  to  the  work  of  Fano 
(Ref  18],  which  was  later  extended  by  Schroeder  and  Aial  (Ref  19]. 
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If  we  replace  the  oscillators  in  the  IPS  definition  by  stable  filters  with  impulse  re- 
sponses of  the  (brni 


h{i)  =  vv(/)^'"^' , 


(82) 


where  w{i)  is  any  real  and  symmetric  function  of  /  with  general  lowpass  characteristics 
and  w(0)  =  1,  then  the  resulting  distribution  Z)(/,  J)  is  given  by 

D,{!,  f)  =  Real  [s\i)is{t)  *  w(/)^'^^']] 

==\[s[l)ls[l)*^,W^''■^'^*  +  /(0[5(/)  *  w( V'^^']] 


(83) 


s{i)s  {x)w{l  -  t)c 


-Jl^Jlr-.)^^  +  -1 


— oo 


5*(/)5(i;)vv(/  -  v)e''^''^'~'^dv. 


Substituting  /,  for  [i  —  r)  in  the  first  integral,  and  —i.  for  (/  —  v)  in  the  second  integral, 
we  ect 


^.fy./)  =  T 


ls[i)s\i  -  /,)  +  s\i)s[i  +  I,)-]  u.-(/,)  e  ^-^'^di 


J 

— oo 


(84) 


=  '^A5..)(^  J) 


where  yir)  =  \v{t  —i). 

That  is.  the  resulting  distribution  is  the  IPS  of  the  windowed  signal,  where  the  win- 
dow is  the  envelope  of  the  filter  impulse  response,  centered  at  the  time  of  analysis.  Since 
\v(j  —  t)  is  real  and  symmetric  around  /,  we  see  from  (  66  )  that  a  smoothing  of  the  IPS 
in  the  frequency  direction  will  result. 

IPS  in  general  can  hence  be  written  as 


IPS{i,f)  =  \ 


is{i)s\i  -  t)  +  s\t)s{t  +  t)]  ^v{0)^v{x)e~^^''■^'dx  , 


(85) 


— oo 


where  \v[r)  is  the  chosen  window. 
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To  illustrate  the  efTectiveness  of  the  windowing  operation,  we  computed  the  IPS  of 
the  test  signal  used  in  Figure  3.  A  Hamming  window  was  used  in  (  85  ),  and  the  results 
are  shown  in  Figure  5.   As  is  observable,  the  ringing  is  avoided. 
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Figure  5. 


IPS  for  a  linear  chirp  -  Hamming  window 


B.     PROPERTIES  AFTER  WINDOWING 

Some  of  the  properties  of  IPS  will  be  affected  by  the  windowing  operation.  It  will 
thus  be  necessary'  to  determine  which  properties  are  modified  and  estabUsh  the  con- 
nections between  the  characteristics  of  the  used  window,  and  the  eflects  it  generates. 

Tracing  back  the  windowing  operation  to  the  kernel  function  ((I>(u,  t))  in  (  20  )  that 
it  implies,  we  will  be  able  to  use  well  estabUshed  results  [Ref  14]  about  these  functions. 

The  resulting  kernel  function  will  be 


0(u>,  t)  =  w(0)vv(t)  cos(7ri)T) , 


(86) 


which  can  be  proved  using  (  19  )  as  follows 


oo        oo       oo 

r    r    r 


IPS{i,  J) 


vv(0)u.(t)  cos(7ri;T)5(/,  +  y  )/(/i  -  y  )e^^''"'e~^^''^"'^'^dodTdt^  .  (87) 


— oo     — oo    — oo 
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Replacing  the  cosine  by  complex  exponentials,  we  have 


oo        oo 


IPSit,  J)  = 


U.-(0)u-fT) 


oo 


5(r,  +  ^  )s\i,  -  ^  )^'''^^'^  ^ T  ^di,e-^'''^'^'^'^dudx 


— oo      — CX5 


— OO 


+ 


oc        oo 


(88) 


%J 

— oo     — oo 


5(/,  +  f  )s\t,  -  -f  )^^'^''<''  -  T  W/j.-^'^'^^'^'^^^^.^T 


— oo 


and,  by  convenient  change  of  variables, 


IPS{i,  f) 


oo        oo 

Vr(0)M(T) 


5(/i)^  (^  —  t)^  '^^^dixe 


'2''J/'j,.^-/2?:(^:+t/) 


i/l)<^T 


— oo     — oo 


— oo 


oc        oc 


+ 


u-(0)u-(t) 


—  >:     — oo 


5*(m)^(a^  +  T)^''"^^^A^e-^''"<^'+^^^u^T 


— oo 


(89) 


Interchanging  the  order  of  integration,  we  obtain 


/?S(/,  f)  =  ^ 


\_s{i)s\i  -  t)  +  5*(/)5(/  +  t)];v(0)vv(T)e"-''^^^Vr 


— oo 


completing  the  proof. 

We  can  now  easily  determine  which  properties  are  maintained  after  the  windowing 
operation,  by  direct  application  of  the  results  in  [Ref  14],  mapping  the  characteristics 
of  the  kernel  function  to  the  properties  of  the  distribution.  The  conditions  on  0(l),  t) 
necessary  for  the  preservation  of  each  of  the  properties  will  also  be  given. 
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1 .  Time  shift 

If  co(/)  =  s{t-  to),  then 

IPSJi,J)  =  IPS,{t-io,J).  (91) 

The  requirement  on  <t>{v,  r)  is: 

<t){v,  t)  constant  for  all  / .  (92) 

This  property  will  thus  be  maintained  for  ever>'  used  window. 

2.  Modulation 

l[coii)  =  s{t)eJ^''^o',  then 

IPSJi,  f)  =  IPS,{t,f-fo) .  (93) 

The  requirement  on  0(0,  t)  is: 

cD(i;,t)  constant  for  all/.  (94) 

This  property  will  thus  be  maintained  for  ever\-  used  window. 

3.  Marginal  in  Time 


oo 


IPS{t,f)df=\s[t)\\  (95) 


— oo 


The  requirement  on  0(i.),  r)  is: 

0(i),  0)  =  1      for  all  V  .  (96) 

Maintained  ifrw(O)  =  1  . 

4.   Marginal  in  Frequency 


oo 


IPS{t,  f)di=\S(J)\\  (97) 


— oo 


The  requirement  on  <I)(u,  t)  is: 

0(0,  t)  =  1      for  all  t  .  (98) 

Not  maintained  for  practical  windows. 


30 


5.   Realness 

IPS{i,  J)  is  real  for  all  /,  /. 
The  requirement  on  0(u,  t)  is: 

O(o,  t)  =  cD*(  -  0,  -  t)  . 


(99) 


(100) 


Maintained  if  the  window  is  real  and  even,  or  at  least  possesses  conjugate  symmetr}' 
around  the  orisin. 


6.   Instantaneous  Frequency 


f-IPS{t,  f)df 


I  s{t)  \' 
The  requirement  on  cL>(u,  r)  is: 


=fi{t)      [instantaneous  frequency) 


(D((>.  0)  =1      and      ^  (D(d,  t)  U.  =  0       for  all  o. 

CT  '~^ 


Maintained  if 

•  ;v(0)  =  1  . 

•  h-'(t)  U  =  0 

7.    Group  Delay 


tIPS{i,f)dt 


J 


SiJ) 


=  ToiJ)      [group  delay) 


The  requirement  on  0(0,  t)  is: 


0(0,  t)  =1      and       ^  cD((;,  t)  ^o  =  ^        ^or  all  r. 


CD 


Not  maintained  for  practical  windows. 
S.   Zero  Power 

5(/o)  =  0  =>  IPS[t,,  j)  =  0 


(101) 


(102) 


(103) 


(104) 


(105) 


This  property  is  clearly  maintained,  due  to  the  multiplicative  nature  of  the  win- 
dowing operation. 


S(^)  =  0  =>  IPS[t,  /o)  =  0 


(106) 


This  property  will  not  be  maintained  since,  from  (  64  ),  a  convolution  along  the 

frequency  axis  will  be  performed. 
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C.      SPECTRAL  RESOLUTION 
{.      Time-invariant  Components 

Though  IPS  is  a  tool  directed  towards  the  analysis  of  signals  with  time-varying 
spectral  contents,  some  insight  into  its  inner  workings  will  be  gained  by  considering  its 
resolution  capabilities  for  thne-invariant  spectral  components.  A  comparison  of  resol- 
ution capabilities  Ibr  IPS  and  WD  will  be  made,  and  contrasted  with  the  standard 
Fourier  Transform. 

As  seen  in  (  6  ).  IPS  is  a  coherent  average  of  two  terms,  one  of  which  uses  only 
past  data,  while  the  other  uses  only  future  information.  When  analyzing  finite  duration 
signals,  the  "past  term"  will  have  its  maximum  resolution  capability  at  the  end  of  the 
data  segment,  since  it  is  at  this  point  that  the  most  past  data  is  available.  Similarly,  the 
"future  term"  will  resolve  better  at  the  beginning  of  the  data  segment.  This  can  be  seen 
in  Figure  6  and  Figure  7,  for  the  past  and  future  terms,  respectively. 


Contour  plot 
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Figure  6.       Sinusoidal  component.  Past  term.  3-D  (a)  and  Contour  (b)  plots 

If  the  two  terms  are  averaged  as  in  (  6  ),  the  resulting  expression  (IPS)  will  keep  the  good 
resolution    at    end-points   that   each   component   provides,   degrading   its   resolution 


A  note  of  caution  should  be  made  here.  These  "past"  and  "future"  terms  should  not  be 
identified  with  the  two  terms  in  (  12  ),  since  each  one  of  the  terms  in  tliis  formula  spans  all  the  data, 
from  —  oo  to  CO. 
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capabilities  towards  the  center  of  the  lime  segment.  This  eflect  can  be  observed  in  the 
contour  plot  of  Tigurc  8. 
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Figure  7.       Sinusoidal  component.  Future  term.  3-D  (a)  and  Contour  (b)  plots 
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Figure  8.       IPS  for  a  sinusoidal  component.  3-D  (a)  and  Contour  (b)  plots 
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This  is  one  of  the  characteristics  of  the  unwindowed  IPS.  It  has  better  resolution  at 
end-points  than  it  does  in  the  middle  of  the  data  segment.  We  should  point  out  that  the 
WD  does  the  opposite,  presenting  a  thinner  main-lobe  at  the  center  of  the  time  interval, 
and  loosing  resolution  capabilities  towards  the  end-points  (Figure  9). 


Contour    pLoL 


en 
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Figure  9.       WD  for  a  sinusoidal  component.  3-D  (a)  and  Contour  (b)  plots 

A  more  formal  discussion  can  be  made  as  follows:  Let  us  denote  by  CW,J[t)  the 
function  oft  whose  Fourier  Transform  is  WD{to,  J)  and,  similarly,  denote  by  CI,J,r)  the 
function  oft  whose  Fourier  Transform  is  IPS{Iq,  f).  Hence, 


c/f;(T)  =  /(/o-YW/o  +  Y) 


1 


(107) 


C/,^(t)  =  —  [5(/o)5   {!q  -t)  +  S   {lQ)s{fQ  +  t)] 


where  s{i)  is  the  stationary'  signal  of  interest.  That  is,  both  IPS  and  WD  create,  from  the 
data,  a  new  "signal",  whose  Fourier  Transform  is  the  wanted  spectrum  (Figure  10).  It 
is  thus  the  effective  duration  of  this  new  signal  that  determines  the  resolution  obtainable 
in  the  frequency  domain. 
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Figure  10.       IPS  and  WD  :  conceptual  diagram 

Let  us  now  assume  that  s(t)  is  only  known  witliin  a  finite  segment  of  Icngtli  L, 
and  compare  the  tlieorctical  resolution  of  IPS  and  WD  at  both  the  end-points  and  the 
middle  of  the  segment. 
a.     End-points 

Since  botli  end-points  are  treated  similarly,  we  will  only  consider  the  start- 
ing point  of  the  known  signal,  which  we  will  denote  by  Iq.  Let  us  consider  that  a 
rectangular  window  of  length  L  was  applied  to  the  signal  (  Figure  11).  It  is  easily  seen 
from  (  107  )  that 


T^O 


=  0 


(108) 
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and  thus,  \VD,jJ)  (Fourier  transform  of  a  delta  function)  is  a  constant.  On  the  other 
hand,  the  windowing  of  the  data  implies  that  a  window  (  vy,(/)  )  of  length  2L  is  applied 
to  Cyr)  (Figure   12). 


Figure  11.       Window  applied  to  the  signal 
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Figure  12.       \Vindo>v  applied  to  CI,J^x) 
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When  we  now  take  tlic  Fourier  transform  of  C/Jt)  ,  we  will  be  taking  the  transform  of 
a  signal  twice  as  long  as  the  original  data.  This  creates  a  main-lobe  twice  as  narrow  as 
the  main-lobe  that  would  be  obtained  if  we  were  to  Fourier  transform  the  data.  How- 
ever, we  know  from  (  107  )  that  C/Jt)  =  Cr,^{-T).  This  then  implies  that,  from  tjie 
obtained  segment  of  length  2L,  only  half  contains  information.  As  a  result,  at  end- 
points,  IPS  presents  a  main-lobe  twice  as  narrow  as  the  one  obtainable  by  Fourier 
transforming  the  data,  but  without  an  equivalent  improvement  in  effective  resolution. 
b.     Middle-point 

We  now  consider  that  the  window  is  symmetrically  placed  around  /g,  as  in 
Figure  13. 
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Figure  13.       Symmetric  niiidow  applied  to  the  signal 

Again,  it  is  readily  seen  from  (  107  )  that  windowing  the  data  with  a  window  of  length 
L  centered  at  i^  imphes  the  windowing  of  CIF,^(t)  with  a  window  of  length  2L.  The  effect 
on  CI,^{r},  however,  is  to  window  it  with  a  window  of  length  L.  This  discrepancy  implies 
that,  in  the  middle  of  the  interval,  the  WD  has  twice  the  effective  resolution  of  IPS.  Let 
us  note  that,  as  discussed  before,  for  both  OFJt)  and  CI,Ji^r),  only  half  of  the  length 
carries  information.  Hence,  their  Fourier-transform  will  have  its  main-lobe  artificially 
narrowed  by  a  factor  of  two,  without  improvement  in  effective  resolution.  In  summary, 
the  WD  has  at  the  middle  of  the  time  interval  the  same  eOective  resolution  that  IPS  has 
at  end-points,  and  hence  the  same  effective  resolution  obtainable  by  Fourier- 
transforming  the  data. 
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This  fact  of  IPS  having  twice  as  poor  resolution  in  the  middle  of  the  time 
interval  than  at  end-points  will  he  important,  since  the  windowed  implementation 
normally  uses  the  window  centered  at  the  time  of  interest.  This  then  means  that  these 
implementations  are  using,  at  each  time,  the  worst  estimate  (in  terms  of  resolution)  that 
IPS  will  provide.  The  model  of  a  typical  windowed  implementation  of  IPS  can  be  seen 
in  Fieure  14. 


Figure  14. 


Typical  implementation  of  IPS 


The  overall  result  is  that,  for  these  implementations,  IPS  has,  by  a  factor  of  two,  lower 
effective  resolution  for  stationary  components  when  compared  with  the  WD. 
2.      Time-varj  ing  Components 

One  of  the  most  noticeable  characteristics  of  IPS  is  the  fact  that  the  width  of  its 
main-lobe  depends  on  the  dynamics  of  the  signal.  A  typical  example  can  be  seen  in 
Figure  15,  where  a  quadratically  chirped  FM  signal  is  used. 
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Figure   15.        IPS  for  FM  signal  -  Quadratic  chirp,  3-D  (a)  and  Contour  (b)  plots 

This  test  signal  lias  its  instantaneous  frequency  increasing  quadratically  with  time.  '1  he 
width  of  the  niain-lobc  is  directly  proportional  to  the  slew  rate.  In  Figure  16,  a  cross 
section  along  the  frequency  axis  of  IPS  for  the  same  test  signal  is  presented.  The  past 
and  future  terms  are  also  indicated. 


CD 


Q_ 
01 


IPS 

Future 

Term 

Post 

Term 

FREQUENCY  RXIS 


Figure   16.       IPS  lor  FM  signal  -  Quadratic  chirp.  Cross-section 
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As  is  seen,  the  width  of  the  main-lobe  of  IPS  results  from  the  imperfect  alignment  of  the 
two  terms:  the  past  term  seems  to  lag  the  true  position  of  the  instantaneous  frequency, 
showing  a  tendency  to  be  positioned  in  locations  previously  occupied  by  the  signal,  while 
the  future  term  leads  the  true  location  of  the  instantaneous  frequency.  This  lag  and  lead 
effect  become  more  severe  for  faster  dynamics,  hence  broadening  the  main-lobe  of  IPS. 
In  an  attempt  to  explain  this  feature,  one  might  be  tempted  to  argue  that  IPS 
lacks  symmetry  in  its  definition.  Despite  being  clearly  insufficient,  as  we  will  see,  it  will 
be  instructive  to  pursue  this  type  of  reasoning  a  little  further,  since  some  additional  re- 
sults will  emerge.  When  computing  CI,J[r),  we  are,  for  each  t  (without  loss  of  generality 
assume  t  >  0  ),  coherently  averaging  two  terms.  One  results  from  the  present  and  past 
signal  history  {x'{t)x{t  +  t)  ),  and  the  other  results  from  present  and  future  signal  histor\' 
(.v(/).v*(/  +  t)  )  .  To  obtain  IPS.  we  will  extract  the  frequency  contents  of  the  sequence  so 
obtained.  However,  if  the  signal  has  frequency  dynamics,  each  one  of  the  terms  will 
contribute  differently  to  the  frequency  contents  o[  CI,^{t)  .  For  example,  if  the  instanta- 
neous frequency  of  the  signal  is  increasing  in  time,  the  contribution  of  the  present-future 
term  {x'{[).x{i  +  t)  )  will  have  higher  frequencies  than  the  present-past  term,  since  it  is 
centered  in  a  region  with  higher  average  instantaneous  frequency.  The  temptation  thus 
arises  to  compensate  for  the  difference  in  the  centers  of  the  two  regions.  This  effectively 
corresponds  to  forcing  the  two  terms  in  Figure  16  to  ahgn  "correctly",  thus  eliminating 
the  broadening  of  the  main-lobe.  To  do  so,  it  suffices  to  introduce  compensating  phase 
terms  in  the  process  of  computing  IPS  (see  .  Alternate  n-ay  of  computing  the  JVigner-Ville 
below).  When  these  compensating  terms  are  introduced,  the  distribution  that  results  is 
the  Wigner-Ville  distribution.  We  can  now  interpret  WD  as:  the  WD  is  the  distribution 
that  results  bv  alisnine  the  two  terms  of  IPS. 


•      Alternate  nay  of  computing  the  Wigner-Ville 


From  the  previous  discussion,  an  alternate  way  of  computing  the  WD  can 
be  derived,  which  allows  us  to  directly  apply  the  WD  to  real  signals  sampled  at  the 
conventional  Nyquist  rate,  thus  avoiding  the  need  for  oversampling  that  would  be 
required  if  the  WD  was  to  be  computed  in  the  usual  way  [Refs.  7,  10].  The  two 
step  procedure  is: 

1.   Compute 

ls{[)s\t  -  r)e~^''''  +s\i)s{i  +  t)^'^'']^^'^^'-''^/  (109) 


F{r,v)  =  -^ 


J 


-oo 
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2.  Compute 


^VDit,  f)  = 


F[T,v)e~^^''^"^'^dDdT  (110) 


— oo    — oo 


As  can  easily  be  recognized,  the  need  for  a  2-D  Fourier  transform  places  the 
cost  of  this  procedure  well  above  the  one  needed  if  an  interpolation  scheme  is  used.  This 
consideration  may  disappear  in  emerging  fields  such  as  optical  signal  processing. 

All  our  previous  discussion  has  been  based  on  a  loosely  defmed  lack  of  symme- 
tr\',  intimately  related  with  the  concepts  of  past  and  future.  As  pointed  out  earlier,  this 
explanation  is  clearly  insufficient.  To  illustrate  why,  a  diflerent  definition  of  "Instanta- 
neous Power  Spectrum"  is  made  in  Appendix  E,  thus  creating  a  new  distribution.  The 
important  point  to  be  made  is:  not  only  are  the  concepts  of  past  and  future  inexistent 
in  this  new  definition,  but  it  also  possesses  perfect  symmetry.  However,  the  resulting 
distribution  is  ver>'  similar  to  IPS,  with  the  characteristic  data-dependent  width.  Since 
lack  of  symmetry  cannot  be  argued  for  this  definition,  the  explanation  for  the  broaden- 
ing of  the  main-lobe  must  be  found  somewhere  else. 

In  Appendix  C,  the  issue  of  uncertainty  is  addressed.  One  important  result  to 
extract  from  there  is:  the  maximum  obtainable  frequency  resolution  when  analyzing 
signals  with  unknown  frequency  dynamics  is  -Jdfjdi .  This  result  is  in  agreement  with  the 
fact,  proven  by  Rihaczek  [  Ref  4],  that  signals  with  strong  phase  modulation  have,  at 
each  time,  their  energy  concentrated  within  a  frequency  band  B^  of  size 


/    1     d'<^{i) 

Bj  =  .Hz—-r^,  (HI) 

V  -'^     dr 

where  0(/)  is  the  signal's  phase.  It  thus  explains  the  broadening  of  the  main-lobe  of  IPS 
in  terms  of  the  expression  of  uncertainty.  Since  the  uncertainty  region  grows  for  in- 
creasingly faster  dynamics,  so  does  the  width  of  the  main-lobe.  Put  in  simpler  terms,  "the 
faster  it  moves,  the  harder  it  is  to  locate".  According  to  the  results  of  Appendix  C,  any 
algorithm  should  behave  in  this  manner,  unless  it  has,  or  assumes  (directly  or  indirectly), 
information  concerning  the  dynamics  of  the  signal.  However,  not  assuming  is  the  key 
rule  for  a  robust  algorithm,  one  whose  performance  does  not  depend  on  the  class  of 
sisnals  beins  analvzed. 
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D.     THE  CHOICE  OF  A  WINDOW 

We  know  from  {  66  )  that,  if  we  pre-window  the  signal  with  a  window  real  and 
symmetric  around  the  time  at  which  we  are  evaluating  IPS,  the  resulting  IPS  will  be: 


IPS^,.^p,  f)  = 


IPS,{t,  v)-IPS^{tJ-  D)dv 


—oo 


(112) 


But 


IPS,(t,f)  =  -^ 


lw{i)\v{t  -  t)  +  w{t)w{l  +  T)]e  ^'^""-^Vr 


— oo 


(113) 


and,  since   vv(/  —  t)  =  \\\i  +  t), 


IPS^ki,  J)  = 


w{t)\v{[  -i-  t)(?  ^^'^■^'dj 


— oo 


(114) 


=  vv(f)fr(/)^''"^' 
=  v^o(0)fr,(/), 

where  u;,(t)  =  vi.(/  +  t),  that  is.  w^  is  the  used  window  shifted  to  the  origin. 

Hence,  if  a  real  and  even  sliding  window  is  used,  the  smoothing  function  in  (  112  ) 
is  the  Fourier  Transform  of  the  window  vv„(t)  .  It  is  interesting  to  compare  this  result 
with  the  Pseudo-Wigner  Ville  (frequency  smoothed  Wigner-Ville),  where  the  smoothing 
function  is  not  the  Fourier  Transform  of  the  window  itself,  but  the  Fourier  Transform 
of  the  square  of  the  window  (in  a  rescaled  frequency  axis)  [Ref.  6].  We  thus  see  that  all 
the  knowledge  on  window  functions  for  the  Fourier  Transform  can  be  directly  applied 
to  the  choice  of  a  window  to  be  used  with  IPS. 

Though  the  smoothing  function  in  (  1 12  )  is  always  real  (we  are  considering  only  real 
and  even  windows),  it  may  be  negative.  This  can  contribute  to  the  presence  of  strongly 
negative  values  in  the  resulting  IPS,  especially  if  the  side-lobe  structure  of  the  window 
is  not  controlled.  This  thus  suggests  that  windows  with  good  side-lobe  suppression 
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should  be  preferred.     Spectral  resolution  is  another  important  issue  to  be  considered 
when  choosing  a  window. 

As  discussed  in  the  previous  section  (Section  III-C).  IPS  has  a  main  lobe  whose 
width  is  proportional  to  ^dfjdi  ,  in  accordance  with  the  uncertainty  principles  of  Ap- 
pendix C.  Hence,  degrading  the  apparent  frequency  resolution  of  the  unwindowed  IPS 
does  not  become  a  relevant  issue,  if  working  with  windows  whose  effective  duration  is 
well  above  the  reciprocal  of  the  uncertainty  region.  The  width  of  this  region  should  thus 
be  the  criteria  for  the  choice  of  window  size,  but  would  require  a-priori  knowledge  of  the 
dynamics  of  the  signal.  In  the  absence  of  this  knowledge,  optimality  will  not  be  achieved. 
However,  as  we  will  see.  IPS  is  a  very'  forgiving  tool,  in  the  sense  that  its  performance 
is  not  seriously  affected  even  for  large  deviations  from  the  optimal  window  size. 
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IV.     EXPERIMENTAL  RESULTS 

To  illustrate  the  types  of  behavior  predicted  in  the  previous  sections,  some  test  cases 
will  be  presented.  For  all  the  cases  the  test  signals  are  128  point  sequences,  with  one 
spectrum  for  each  sample.  All  plots  are  on  a  linear  scale.  Some  considerations  relating 
to  the  Pseudo  Wiener- Ville  (PWD)  were  made:  all  test  signals  are  analytic,  thus  freeing 
the  PWD  from  the  somewhat  annoying  interferences  between  positive  and  negative  fre- 
quencies. The  exception  to  this  rule  will  be  Figure  19,  since  its  sole  purpose  is  to 
illustrate  the  eflects  of  the  use  of  real  signals.  Also,  when  noise  is  present  (Figure  24), 
analytic  noise  is  used.  L'nless  otherwise  stated,  all  Figures  use  a  41 -point  Flamming 
window,  and  present  the  distributions  after  smoothing  along  the  time  direction.  Hence, 
most  of  the  spectral  cross-terms  of  the  PWD  due  to  multi-component  signals  do  not 
appear  in  the  plots.  To  illustrate  the  edectiveness  of  time-smoothing,  Figure  17  was 
provided,  diflering  from  Figure  21  only  in  the  fact  that  the  time-smoothing  technique 
was  not  used. 


FREQUENCY  AXIS 
(a) 


FREQUENCY  AXIS 
(b) 


Figure  17.       FSK.  IPS  (a)  and  PWD  (b).  No  time-smoothing. 

A.     REAL  SIGNALS 

In    Figure   18,    IPS    and    PWD    for   a    linear   analytic    FM    chirp   are   presented. 
Comparing  Figure   18  with  Figure   19,  where  the  test  signal  is  a  linearly  chirped  cosine 
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function,  wc  can  appreciate  how  insensitive  IPS  is  to  the  fact  of  tlie  siiinal  being  real. 
The  artifacts  prcscin  in  the  PWD  wiien  analyzing  real  signals  are,  for  all  practical 
purposes,  absent  in  IPS. 
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Figure   18.       Analytic  linear  FM  chirp.  IPS  (a)  and  PWD  (b). 
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Figure   19.       Linearly  chirped  cosine.  IPS  (a)  and  PWD  (b). 
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B.  INSTANTANEOUS  POWER 

Another  visible  ellect  in  Figure  19  is  the  "amplitude  modulation"  of  the  main-lobe 
of  IPS.  IPS  gives  us  the  distribution  in  frequency  of  the  instantaneous  power  of  the  sig- 
nal. It  then  follows  that,  when  computing  the  spectrum  at  a  low  power  sample,  the 
amplitudes  of  the  frequency  components  will  have  to  be  small,  since  their  total  power 
must  add  up  to  the  low  power  of  the  sample.  This  is  a  direct  consequence  of  (  95  )  and 
is  a  requirement  that  has  be  fulfilled  by  any  positive  distribution.  It  is  not  observable  in 
Figure  18,  since  the  used  analytic  signal  has  the  same  power  at  all  samples. 

The  WD,  though  also  obeying  the  Marginal  in  Time  property,  tends  to  absorb  the 
fluctuations  of  the  instantaneous  power  in  its  cross-terms,  and,  in  order  to  do  so,  is 
forced  to  assume  strongly  negative  values. 

C.  END-POINT  RESOLUTION 

Though  the  better  resolution  at  end-points  provided  by  IPS  is  in  general  well  visible 
in  all  the  plots  that  are  presented,  Figure  20  is  probably  one  of  the  most  representatives. 
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Figure  20.       FSK  (127-point  Hamming  window).  IPS  (a)  and  PWD  (b). 
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Considering  it,  \vc  sec  that  not  only  is  the  end-point  resolution  of  ITS  better,  but  also 
the  transition  times  arc  shorter,  allowing  a  more  accurate  definition  of  the  time  ol  spec- 
tral jumps. 

As  stated  in  Section  lll-D,  li'S  is  a  forgiving  tool,  in  the  sense  that  we  can  deviate 
from  the  optimal  length  of  window  without  severely  degrading  its  performance.  To 
illustrate,  Figure  21  is  presented,  which  dilfers  from  figure  20  only  in  that  it  uses  a 
window  with  dilfcrent  effective  duration. 


rREOUENCY  AXIS 
(a) 


FREQUENCY  RXiS 
(b) 


Figure  21.       FSK  (  Sl-poinl  Hamming  niiidow  ).  IPS  (a)  and  F^^  D  (b) 

As  is  seen,  the  time  resolution  is  basically  kept  the  same,  the  only  noticeable  difference 
being  the  width  of~the  window's  main-lobe. 

D.     MULTFCOMPONENT  SIGNALS 

Performance  for  multi-component  signals  is  addressed  in  Figure  22  and  Figure  23, 
where  the  test  signal  has  three  components:  a  linearly  chirped  FM  signal,  a  quadratically 
chirped  FM  signal,  and  a  stationary  component.  As  seen,  IPS  is  free  of  the  cross-terms 
appearing  in  the  PWD.  However,  for  IPS,  each  one  of  the  spectral  components  has  its 
width  aflected  by  the  associated  uncertainty  region  and  is,  hence,  wider  than  the  corre- 
sponding one  in  the  PWD.  Again,  the  instantaneous  power  modulation  is  visible  in  IPS. 
These  effects  are  also  apparent  in  the  provided  contour  plots. 
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Figure  22.       IPS  for  Multi-component  signal.  3-D  (a)  and  Contour  (b)  plots 
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Figure  23.       PWD  for  Multi-component  signal.  3-D  (a)  and  Contour  (b)  plots 
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E.     PERFORMANCE  IN  NOISE 

One  issue  still  not  addressed  is  the  performance  in  noise.  In  Appendix  I),  a  proof  is 
given  that,  for  a  signal  embedded  in  additive  bandlimitcd  uhite  Gaussian  noise,  the 
variance  of  IPS  is,  on  the  average, -''  smaller  than  the  variance  of  WD.  Also,  from  the 
same  Appendix,  we  have  that,  wiien  the  bandwidth  of  the  noise  approaches  infinity,  the 
improvement  in  the  variance  given  by  II'S  approaches  the  limiting  value  of  3  dB. 
Figure  24  is  provided  to  illustrate  this  theoretical  result.  The  SNR  for  this  picture  is 
5  dB,  and  it  is  seen  that  IPS  does,  indeed,  present  a  smoother  estimate  of  the  instanta- 
neous spectrum. 
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Figure  2-1.       I  SK  (  SNR  =  5  dD  ).  IPS  (a)  and  PVVD  (a) 

F.     THE  CAPTURE  EFFECT 

Having  borrowed  the  term  from  communication  theory  (FM  demodulation),  our 
goal  is  to  compare  the  behaviors  of  II\S  and  PWD  when  analyzing  multi-component 
signals,  when  the  components  have  diflcrent  amplitudes.  In  Figure  25,  the  two  compo- 
nents have  amplitudes  dificring  by  a  factor  a  four  (12  dB).  As  is  observable,  both  IPS 
and  PWD  detect  the  weaker  component,  despite  the  fact  of  the  cross-terms  in  the  PWD 
being  already  more  energetic  than  the  weak  signal  itself    In  Figure  26,  the  same  two 


■^  The  variance  of  IPS  is  time  dependent.  When  we  say  that  "...the  variance  is. ..on  the  aver- 
age....", wc  arc  rclcrring  to  the  time  average  of  the  variance. 
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component  signal  is  presented,  now  with  a  factor  of  8  between  amplitudes  (18  dB).  It 
is  obvious  that  the  PWD  lost  the  weak  component.  The  cross-terms,  however,  still  re- 
main. In  this  case,  the  position  of  the  weaker  term  could  be  inferred  from  the  cross- 
terms,  but  such  a  technique  would  be  impossible  to  apply  to  more  complex  signals. 


PREOUENCY  AXIS 
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Figure  25.       Po>ver  ratio  =  12  dB.  IPS  (a)  and  PWD  (b) 


FREQUENCY  AXIS 
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Figure  26.       Power  ratio  =   18  dl3.  IPS  (a)  and  PWD  (b) 
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V.     PARAMETRIC  IPS 

A.  INTRODUCTION 

We  have,  from  (  12  ),  that  the  IPS  of  a  signal  is,  at  each  time,  the  Fourier  transform 
of  a  function  of  the  signal.  This  imphes  that  the  wanted  information  hes  in  the  spectral 
content  of  this  function.  The  process  of  obtaining  the  IPS  can  therefore  be  reduced  to 
the  problem  of  performing  conventional  spectral  estimation  on  a  function  G{{,  t)  which 
is  derived  from  the  signal  as  in  (  II  ).  The  use  of  parametric  methods  should  thus  be 
considered,  in  an  attempt  to  improve  the  resolution  capabilities  of  IPS.  These  consider- 
ations are  also  applicable  to  the  WD,  and  a  parametric  version  of  this  distribution  has 
been  considered  [Ref  20]. 

B.  AR  MODELING  -  CONSIDERATIONS 

Though  obtained  as  a  bilinear  transformation  of  the  data,  the  function  in  (  11  )  can 
not  be  considered  a  true  autocorrelation  function,  since  it  is  not  constrained  in  any  way 
other  than  possessing  conjugate  symmetr\'  around  the  origin.  It  can  not  be  guaranteed 
to  give  rise  to  a  positive  definite  autocorrelation  matrix,  and  hence  can  not  be  used  to 
directly  solve  the  normal  equations.  We  can  however  estimate  its  autocorrelation  func- 
tion and  fit  the  model  to  this  autocorrelation.  In  a  sense, '^  we  will  be  fitting  the  model 
in  the  fourth  moment  of  the  data.  In  all  the  results  presented,  the  AR  parameters  are 
obtained  using  the  Modified  Covariancc  .Method.  This  method  has  been  found  to  pos- 
sess good  resolution  properties,  while  alleviating  some  of  the  problems  associated  with 
the  Maximum  Entropy  .Method  (.MEM)  [Ref  21).   All  plots  are  on  a  logarithmic  scale. 

C.  EXPERIMENTAL  RESULTS 

Choosing  the  order  of  the  model  is  a  sensitive  issue.  From  Figure  16.  we  know  that 
the  main  lobe  of  IPS  consists  of  two  terms,  which  do  not  align  exactly.  We  might  thus 
expect  that,  if  the  correct  order  is  exceeded,  the  extra  poles  will  tend  to  resolve  these  two 
terms,  and  each  spectral  component  will  have  two  associated  peaks  in  the  spectrum.  In 
Figure  27,  where  the  test  signal  is  again  a  single  component  quadratically  chirped  F.M 
signal,  this  elTect  is  well  visible.  In  the  left,  a  one-pole  model  was  used  (analytic  test 
signal).  With  the  addition  of  one  extra  pole,  we  can  see  that  the  two  terms  were  resolved 


Lacking  both  temporal  and  statistical  averaging.  (11)  can  not  be  considered  to  represent  the 

second  moment  of  the  data. 
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(with  some  bias),  and  the  undcsired  eflect  is  now  present.    In  Figure  28,  a  fourth  order 
model  was  used  to  fit  a  three  component  signal. 
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Figure  27.       Parametric  IPS.  First  (a)  and  second  (b)  order  models 

We  can  observe  that,  since  the  main  lobes  were  not  over-resolved,  the  signal  components 
are  well  defined,  and  can  be  located  with  greater  precision  than  what  would  have  been 
possible  with  the  non-parametric  IPS. 


Figure  28.       Tliree  component  signal.  Parametric  IPS.  Fourtli-order  model. 
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VI.     CONCLUSIONS  AND  RECOMMENDATIONS 

In  this  thesis,  an  ahernative  tool  for  the  spectral  analysis  of  signals  with  time- 
vaning  spectral  content  is  presented  (IPS).  Having  its  roots  in  the  definition  proposed 
by  Page  of  "Instantaneous  Power  Spectrum",  IPS  has  a  defining  expression  ver\'  similar 
to  the  one  defining  the  Wigner-Ville  distribution  (WD).  Their  performance  is,  however, 
considerably  difierent.  The  WD  approach,  though  extensively  used  in  the  last  decades, 
suffers  fi-om  a  number  of  short-comings.  The  foldover  frequency  of  its  discrete  fre- 
quency is  located  at  nj2  ,  requiring  the  use  of  analytic  signals  or  the  use  of  an  interpo- 
lation scheme  as  an  alternative  to  sampling  at  twice  the  Nyquist  rate.  Interference  when 
used  with  real  signals,  cross-terms  when  multi-component  signals  are  present,  preference 
of  linear  dynamics  and  loss  of  resolution  at  the  extremes  of  the  analysis  segment  are 
other  problems  associated  with  the  WD  [Ref  8].  The  distribution  presented  here  is 
shown  to  have  some  advantages  over  the  WD.  These  are:  evenness  of  performance  for 
different  dynamics  of  the  signal,  direct  applicability  to  real  signals  sampled  at  the 
Nyquist  rate,  reliability  when  analyzing  multi-component  signals,  and  preservation  of 
spectral  resolution  at  end-points.  Also,  IPS  is  shown  to  behave  better  in  the  presence 
of  additive  white  Gaussian  noise  (AWGN).  As  a  disadvantage,  IPS  does  not  achieve  the 
frequency  resolution  provided  by  the  WD.  Due  to  its  robustness,  IPS  is  especially  well 
fitted  to  be  used  as  a  front-end  tool  in  non-stationary  spectral  analysis.  The  usefulness 
of  AR  modeling  as  applied  to  IPS  is  investigated.  It  is  found  to  be  a  technique  ver\' 
sensitive  to  the  chosen  order  model,  but  providing  some  advantage  in  defining  the  lo- 
cation of  each  spectral  component.  The  use  of  other  parametric  techniques  should  be 
investigated.  In  Figure  29.  we  can  see  the  FSK  signal  of  Figure  20,  this  time  processed 
with  the  following  expression: 


D{t,J)  =  4^ 


lsii)s\t  -  r)  -  s\t)s{t  +  t)]^"-^^^-^Vt  .  (115) 


— oo 


That  is,  Figure  29  is  the  imaginary  part  of  the  Rihaczek  distribution  of  the  test  signal. 
As  is  seen,  it  shows  some  natural  ability  to  detect  fast  transitions  in  the  spectrum,  and 


53 


seems  to  be  a  promising  area  for  future  work.   Also,  linear  combinations  between  IPS 
and  this  imaginary  part  of  the  Rihaczek  distribution  should  be  investigated. 
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Figure  29.       Imaginary  part  of  Rihaczek's  Distribution  -  FSK 

Signal  synthesis  from  IPS  may  also  be  a  promising  field  of  research  for  speech  processing 
applications. 
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APPENDIX  A.     IPS  FOR  FINITE  DURATION  DISCRETE  SIGNALS 

In  this  appendix,  a  formal  derivation  of  the  expression  defining  IPS  for  discrete 
signals  {  16  )  is  presented,  as  an  alternative  to  the  direct  discretization  of  (  12  ). 
The  discrete  version  of  (3)  is.  assuming  that  the  sequence  starts  ixt  n  =  0: 


S-{n,e)  =  AT)    s(n}e-^^"  (116) 

n=0 


Hence, 


where 


and 


Similarly 


S~in,  8)  1^  =  AT' 


y    s{n)c  -^  '/    s  [n)e'  " 


k=-i 


n=i 


(117) 


c-=  \    s{n)s*{n-k)       k>0  (US) 


c:  =  {c:S-  (119) 


,V-1 

S^{n,d)  =  ATy  s{n)e-'^'  (120) 
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^  ,v-i 


|sV^)P  =  Ar^)  5(«).-^''"  ^ /(«y^'' 


n=t 


(-V-1-0 

z 

*=-(-V-i-0 


(121) 


=  Ar^    ^     c,+  .-^-^^ 


where 


and 


c+=    >     5(a7)5*(«-A)       k>0  (122) 


c,+  =  (0*.  (123) 


Now,  to  get  the  digital  equivalents  or(2),(4)  we  must  approximate  the  derivatives  of 
(   117  )  and  (  120  ),  which  we  will  do  using  first  order  diflerences: 


p~{}u  6)^  ^T[  I  S~{ih  0)\^  -\  S~{t  -  i,e)f] 


k=-i 


(124) 


where 

zj.  =  s{n)s\n  -  k)      k>0  (125) 

and 

^I  =  iCS  (126) 

In  the  same  way, 

pV  e)^AT[\  s>.  e)\^  +  \  s^{i  +  i,  e)  \^] 

(■V-1-0 
AT        V  -.     -JOK  (^27) 

k=-{\-]-!) 
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where 

4  =  s\n)s{n  +  k)      k>0  (12S) 

and 

4  =  (0*-  (129) 

Assuming  the  signal  to  be  zero  outside  the  known  samples,  then  according  to  Levin's 
definition  (6): 

AT     V^  .  ^^^''^ 

=  -^    /     [s{n)s\n  -  k)  +  s\n)s{n  +  k)-]e-^^^ 
k=—oo 

which  is  the  wanted  result. 
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APPENDIX  B.     INSTANTANEOUS  FREQUENCY  AND  IPS 

The  proof  of  property  6  (  55  )  can  be  made  as  follows: 


r 


/■IPS{[,  M 


J 

— oo 


^(0  r 


(131) 


oo  oo 


5(0 


f— 

^  2 


\_s(i)s\t  -  t)  +  ^\i)s{t  +  x)-]e~^^''^'dxdf 


— oo  — oo 

oo 


2  1  5(0 


2  1  5(0 


r 


[5(05    {t-x)  +  S   {1)S{1  +  t)] 


fe-J^'-^'dfdT 


— oo 


%J 

— oo 


[5(/)5*(r  -  t)  +  s\t)s{i  +  t)]  -r^  -f-  [^ Wl^^ 

2?:    at 


1  /     J 


2  I  5(0 


2  \    :>-T 


5(05  (/  -  T)d  {T)dr  + 


(132) 


+ 


1  /    J 


2  I  5(0 
J 


5  (/)5(/  +  T)d  {rjdz 


— OO 


y 


-5(05(0- ^ -J- 5   (05(0 

4:r  1  5(0  I'  4;r  |  5(0  |' 


^         -[5(05*'(0-5*(05'(0] 

-[-y2Imag[5*(05'(0]] 


An  I  5(0 


An  1  5(0 
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Therefore 


fips{t,  M 


5(/) 


Imag[5  {i)s  (/)] 


(133) 


For  real  signals. this  cenier  of  mass  is  zero,  as  expected,  since  the  spectrum  is  conju- 
gate symmetric  about  DC. 
For  analytic  signals: 


s{t)  =  m{i)e''^^'^ 


(134) 


and 


f-iPS{i,  M 


5(0 


Imag[5  (/)5(0] 

271  I  5(0  l' 

Imag[mfOe--^''^'^^[/n'(O^'''^V70'(Om(Oe^''^^^^]] 
2nm^{i) 

Imag[i7;(0m  (/)  +./0  (0'^''(0] 
27im'^[i) 


(135) 


and.  fmallv 


m=i^im)^ 


(136) 


which  is  the  definition  of  instantaneous  frequency  for  analytic  signals. 


This  defuiilion  of  Instantaneous  Frequency  was  introduced  by  Ville  [Ref.  2]  in  1952. 
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Hence,  we  have  that,  for  analytic  signals 


and  our  proof  is  complete. 


f-IPS{t,  f)df 


S{1) 


-m 


(137) 
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APPENDIX  C.     NON-STATIONARY  SIGNALS  AND  THE  PRINCIPLE  OF 

UNCERTAINTY 

When  analyzing  signals  with  time-varying  spectral  contents,  there  is  a  maximum 
obtainable  frequency  resolution,  determined  by  the  dynamics  of  the  signal.  This  is  in 
contrast  with  the  stationary  case,  where  the  only  constraint  on  the  obtainable  resolution 
is  the  observation  time.  The  presentation  and  proof  of  this  effect  is  the  goal  of  this 
chapter. 

A.      THE  PHYSICAL  CONCEPT 

Before  presenting  any  mathematical  results,  we  will  try  to  approach  the  subject  with 
a  brief  discussion  of  the  involved  ideas. 

Let  us  consider  our  abihty  to  perform  spectral  estimation  on  a  stationary'  signal. 
For  that  purpose,  we  will  assume  that  the  incoming  signal  has  been  segmented  into 
contiguous  intervals  of  A/  seconds.  Our  final  goal  is  to  determine  the  spectrum  occupied 
by  the  signal  at  roughly  the  moment  of  analysis,  that  is.  within  the  last  received  segment. 

If  we  are  given  only  this  last  segment,  the  obtainable  frequency  resolution  is  con- 
strained by  [Rcl'.  1] 

A/>    j^  (138) 

But  if  now.  for  some  reason,  the  next  to  the  last  segment  of  data  becomes  available. 
we  will  be  able  to  improve  our  analysis  of  the  last  segment  up  to  a  resolution  of 

We  can  continue  the  process  and,  if  more  segments  of  data  become  available,  im- 
prove our  spectral  description  of  the  signal  at  the  present  time.   The  rule  in  this  case  is: 


•  The  stationarity  of  the  signal  allows  us  to  improve  the  analysis  of  the  present  time 
by  using  data  collected  no  matter  how  distant  in  the  past.  Hence,  if  infmite  obser- 
vation time  is  allowed,  mlinite  spectral  resolution  can  be  achieved. 
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Let  us  now  consider  the  same  scenario,  but  with  a  non-stationar\'  incoming  signal, 
v/hich.  we  will  assume,  can  be  considered  stationary-  within  the  segment.  If  only  the  last 
segment  is  available,  our  frequency  resolution  is  again  constrained  by 


A/> 


1 


2-Ar 


(140) 


and  the  only  way  to  improve  it  is  to  consider  more  segments.  But  a  qualitative  difference 
appears  at  this  point.  Segments  of  the  "distant  past"  are  rendered  useless  by  the  dy- 
namics of  the  signal  unless  we  know  or  assume  its  behavior.  Considering  those  segments 
will  only  degrade  our  spectral  analysis.  The  usable  data  is  thus  restricted  to  the  most 
"recent  past",  which  implies  an  upper  bound  on  the  obtainable  resolution. 


•  For  a  non-stationary  signal  with  unknown  dynamics,  there  is  a  maximum  obtaina- 
ble frequency  resolution,  which  depends  on  how  distant  the  unusable  "distant 
past"  is.  It  is  an  absolute  maximum,  and  depends  only  on  the  dynamics  of  the 
sisnal. 


Proving  this  statement  is  the  purpose  of  this  Appendix. 

B.     MATHEMATICAL  FORMULATION 

Let  us  briefly  discuss  Gabor's  result.  Gabor  [Ref  1]  defined  the  "effective  duration" 
A/  and  the  "effective  frequency  width"  A/ of  a  signal  4/(1)  by  the  following  equations 


Ai=^^2Ai-i 


-.2 


(141) 


Af=^2r:{f-jy 


where  [.]  stands  for  the  average,  and  the  «'''  moments  are  defined,  ommiting  the  argu- 
ment of  >p{i).  as 


<A*/V# 


in = 


r 


\p  \pd[ 


— oo 


C/')  = 


2nj 


at 


i//  xj/dt 


— oo 


(142) 
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resulting  in 


A/'A/> 


(143) 


This  only  states  that  any  signal  whose  eflective  duration  is  At  (  sec  Figure  30  ),  will 


have  an  eflective  bandwidth  A/(  see  Figure  31  )  greater  than 
tainty  is  not  addressed  in  this  result. 


1 


2A/ 


;  the  issue  of  uncer- 


1 

\/yCt) 

^ 

At 

Figure  30.       Effective  duration  of  s(t) 


Figure  31.        Effective  bandwidth  of  s(t) 

In  fact,  its  relevance  as  an  uncertainty  principle  only  comes  into  play  when  we  realize 
how  classical  Spectral  Estimation  is  achieved. 

Classical  methods  are  based  on  the  "measure  of  similarity"  between  the  data  and 
more  elementary  functions  of  time  which  have  strong  affniities  with  the  physical  concept 


63 


of  frequency.  The  best  example  is  the  Fourier  Transform,  which  conceptually  is  a  bank 
of  correlators,  tr\'ing  to  determine  the  "measure  of  similarity"  between  the  data  and 
sinusoidal  functions,  the  true  representation  of  our  concept  of  frequency.  Other  typical 
examples  are  in  the  works  of  Gabor  [Ref  1]  and  Priestley  [Ref  22],  using  elementar}' 
functions  on  the  t  —f  plane,  where  these  functions  have  what  Priestley  calls  an 
"oscillatorv'  form"  needed  to  establish  the  correspondence  to  our  notion  of  frequency. 

Now,  if  the  data  is  of  length  T,  the  elementar\^  functions  can  have  at  best  the  same 
duration,  which  implies  that  their  frequency  width  is  greater  than  -rf.  Since  the  ob- 
tainable frequency  resolution  is  at  best  of  the  size  of  the  frequency  width  of  the  ele- 
mentar\'  functions,  the  principle  of  uncertainty  in  Spectral  Estimation  follows. 

Gabor's  result  is  only  concerned  with  stationary  signals,  or  at  least  assumed  to  be 
stationary',  since  A/ is  the  effective  width  of  the  Fourier  Transform  of  the  signal.  How- 
ever, some  imphcations  can  also  be  derived  for  the  non-stationary  case. 

Let  us  assume  we  have  an  arbitrary"  stationary  signal  \li{i),  which  we  will  consider  to 
be  analytic  to  simplify  the  formulation,  with  effective  duration  A/,  effective  bandwidth 
A/ and  constant  instantaneous  frequency.    Formally, 

1      d 


^_         \_dir%\_\ljU)']^  =  constant  (144) 

Let  us  now  force  the  signal  to  linearly  chirp  in  frequency,  by  multiplying  it  with  the 
phasor  e''''"'^  and  defining  the  modified  signal  i//,w(^)  by 

where  A  is  a  positive  real  number.  Due  to  the  exponential,  ^J^t)  is  now  a  non-stationar\' 
signal  chirping  linearly  in  frequency  with  a  slew  rate  of 

±-k.  (146) 

The  following  results  can  be  established  for  ^m{^)'- 

•      Effeciixe  lime  duration 

The  effective  time  duration  will  be 


^t,,  =  ^2n{l^--l^Y  (147) 
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where 


Hence, 


^M^M^dt 


^M^'M^l 


— oo 


oo 


ij/  ip  t  dt 


4/  ij/dt 


=  /• 


— oo 


^M^I^M^    ^^ 


\p  ij/ 1  dt 


2  \         — oo 


(',y = 


— oo 


= ('^)  • 


^h]^hidt 


— oo 


(//  t/'^/' 


— oo 


(Afv,)'  =  271  [  (/')  -  (F)']  =  {Atf 


•      Effective  bandwidth 

The  efTective  bandwidth  will  be 


where 


(148) 


(149) 


(150) 


AA,  =  ^2n(J\r  -]\,f  =^  27r  [  (f,,f  -  (f\ff  ]  •  (151) 


^*M-^^M^^ 


f\i  —    -) 


1        ±-_ 


-^7         z' 


^M^M^^ 


(152) 


Using  the  identity 


.,     =f[.W.'-1  =  f^-Vi2,A.,^-=, 


(153) 
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we  have  that 


"A 


M 


Im- 


1        *^ 

1  — oc 


^  /^"^''  +  /^2;rA7^^-'^^'' 


^/ 


dt 


2nj 


\p  ij/dt 


— oo 


oo  oo 

I-  '''I'  ^ 

i,  —d, 


1     ^ 

i         — oo 


i/'  \l/kidi 


InJ 


+ 


— oo 


\p  \j/dt 


— oo 


l//   \pdl 


oc  oo 


1     *^ 

1  —  rxD 


(//  l/'^'^'r 


2^7 


+  /c 


—  DO 


i//  \\)dt 


l//  l//<^/' 


(154) 


Hence,  all  we  are  missing  to  completely  determine  A/^,  is  the  expression  for  [j\) 
where  [Ref.  1] 


(^v5) 


d,        d,     '" 

1     ^- 

-oo 

{Inf 

OC 

%. 

-oo 

(155) 


Noting  that 

-^-^  =  —  -^-^2"^^^  — +,2./cr^^^  +  (2.A')/^^   .    (156) 
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we  have 


^■^-^ 


[It: 


[jj  ij/di 


— oo 


(157) 


(//  I  \pdt 


=  if)  +  k 


2    ^>r- 


+ 


Jk 


ibi — ; — dt 
di 


Jk    i 


^i  —  dt 


An 


ip  ipdl 


—  30 


Zn  ^ 


}p  \pdl 


— oo 


l//  ijjdl 


— oo 


Con<;idering  now  the  following  identities 
dip*        i//tA*    d 


d 


*  d\J/         \Ij  i//     //  *   d 


we  have 


(158) 


(A/./  — [arg(i//)](y/ 


t/;})  =  (/i  +  A^(rVi 


A'        — >3 


i//  i/z^r 


— oo 


(159) 


and,  using  (  144  )  we  get 


{f.i)  =  if)  +  ky)  +  2kif, 


(160) 
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Since  [Ref.  2] 


we  have,  by  (  144  )  that 


ij/fiil/dt 


/  = 


1     ^ 

1        — oo 


\\)  \lidi 


— oo 


(161) 


From  where, 


i//iA  dt 


f=fi 


=  fi 


\p  \pdt 


— oo 


(/^,)  =  if)+ky)  +  2kjf. 

Hence,  from  (  151  )  and  (  154  ). 

A/v/  =  2;r[(/-)  +  kV)  +  2Ai/]  -  2n[(Jf  +  k\lf  +  IklJ] 
=  2n[Jf)-iJf]^k'[i/)-ilf] 


=  {Af)'  +  k\Aif  . 


(162) 


(163) 


(164) 


Therefore,  using  (  150  )  we  can  write  the  effective  time  and  frequency  widths  of  the 
modified  signal  \1/m{^)  in  terms  of  the  effective  widths  of  the  original  stationary  signal 
ij/it)  as 


{Af^f  =  {AJ)'  +  k\Aif 
(Ai,,)'  =  [Atf  . 


(165) 


Using  (  143  )  and  (  165  ),  we  get  the  result 
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{Af,fAr,^)'  =  {Aif[{Af)'  +  k'{Aif] 


2/  *  .x4 


=  (A/A/)-  +  k\Ai 
>-^  +  k\Aif, 


(166) 


from  where 


4{Aif 


(167) 


Hence,  if  we  minimize  the  right  hand-side  of  (  167  ),  we  will  find  a  lower  bound  on  the 
efTective  bandwidth  of  the  chirped  signal.  Proceeding  with  the  minimization,  we  have 


<5A/  1_    4(Ao 


^         +k\A!f 


=  0 


4(A/) 


+  2A-  (A/)  =  0 


>k^A[ 


4{Aiy 


4k' 


=^At  = 


V   2k 


Hence 


{Af,,f  >  4 + 4 


(169) 


AA/  ^  %■ -^ 


and  we  can  read  the  obtained  results  as  stating: 


Lemma  J.  A  linear  chirp  will  always  occupy  an  effective  bandwidth  greater  than 


y:df,jdi  ,  independently  of  its  efTective  duration. 
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Lemma  2.  A  linear  chirp  will  occupy  the  minimum  efTective  bandwidth  o^  ^jd/Jdi , 
only  if  its  efiective  duration  is 

At  = 


These  results  will  have  obvious  impHcations  when  classical  spectral  analysis  of  the 
chirped  signal  is  attempted;  no  matter  what  the  duration  of  the  elementary  functions  is, 
the  signal  will  always  appear  to  occupy  more  than  -/d/Jdi ,  and  this  will  thus  be  the  best 
obtainable  resolution.   This  completes  our  proof 

Lemma  i.The  maximum  obtainable  frequency  resolution  when  analyzing  signals 
with  unknown  frequency  dynamics  is  ^dfjdi  . 

Though  the  previous  discussion  was  restricted  to  the  case  of  a  linear  chirp,  the 
underlying  ideas  remain  valid  for  laws  of  variation  other  than  linear;  the  treated  case  can 
be  considered  as  a  first  order  approximation  for  arbitrary  non-stationary  signals. 


70 


APPENDIX  D.     PERFORMANCE  IN  NOISE 

As  we  will  see  later  in  this  Appendix,  an  absolute  evaluation  of  the  performance  of 
IPS  in  noise  is  not  possible,  because  the  second  moment  is  data  dependent.  We  will  thus, 
instead,  compare  the  noise  performance  of  IPS  with  the  noise  performance  of  the 
Wigner-Ville  distribution.  Let  us  consider  the  inverse  Fourier  transforms  of  IPS  and 
WD. 


C/,(r,  t)  =  ^  ls{t)s\t  -t)  +  s\i)s{t  +  t)] 

versus  (172) 

aV,{t,r)^s{i+  j-)s\t-j-), 

Given  a  signal  z{i}  that  is  corrupted  by  WAGN,  where  the  noise  is  analytic  and 
bandlimited  to  the  maximum  frequency  (H")  of  the  signal,  that  is: 

v(/)  =  r(/) +  «(/),  (173) 

where 

r.Vn        0  <  u-  <  IV 
5„(vv)  =  <  (174) 

0  otherwise 


1 


A.    MEAN 
1.      WD 


Cn\{[,  t)  =  v{[  +  j-  )v\r  -j-)  =  [z{[  +  j-)  +n{i  +  j-  )]-[z\i  -j-)  +  n\i  -f)] 

=  z{t  +  f  )z\i  -^)  +  z{t  +  f  )n\i  -  j- )  (176) 
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Hence, 


E{CU\i[,  T)}  =  aV,iL  t)  +  zit  +  f  )E[n\i  -j-)] 

+  z\i-j-)E[n{t  +  j-)n\t-j-)] 

E{Cn\{i,T)]  =  aV,{t,r)  +  R„{r). 


(177) 


2.      IPS 


Clyit,  T)  =  ^  [  v{i)v\t  -  t)  +  v\t)v{l  +  t)  ] 

=  ^  [z(/)  +nit)yi2\t  -  t)  +  «*(/  -  t)]  (179) 

+  -^\_z\i)+n\imz{l-r)  +  n{i-T)-]. 


E{CI,{t,  t)  }  =  -f  z(/)r*(/  -  t)  +  ^  E{n{t)n\t  -  t)} 


2 
'2 


+  -^z*(02(r  +  t)  +  4"  4^*(0«(^  +  t)}  (180) 


=  CL{i,t)  +  \r,{t)  +  ^R,{t) 


Hence, 


£(C/,(/,t)}  =  C/,(m)  +  /?„(t).  (181) 


Therefore,  IPS  and  WD  have,  as  far  as  their  first  moment  is  concerned,  identical  be- 
havior in  the  presence  of  noise. 

B.     VARIANCE 
1.     WD 


I  CW,{i,  t)  I'  =  v(/  +  ^  )v*(^  -  -^  )v\t  +  ^  )v(/  -  -^  ) 

=  [z(:  +  ^)  +  «(/  +  -^)][z*(/-y)  +  «V-y)] 


.[/(,  +  ^ )  +  ;,•(,  +  -I- )][,(,  _  ^  )  +  „(,_  ^ )] 
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+  n{i  +  j-  )z\l  -j-)  +  n{i  +  j-  )n\t  -  ^  )] 

•[2  V  +  j-  )z{t  -j-)  +  z\i  +  j-  K/  -  Y ) 

Using  now  the  expectation  operator,  we  have: 

E{\ClV,{i,T)\']  =  \z{!  +  j-)\'\z{i-j-)\' 

+  z{t  +  ^)z\t-j-)E[n\!  +  j-)] 
+  \z{i  +  j-)\'E[\n{t-j-)\'] 
+  z{i  +  ^  )z{!  -  f  )E\n\i  +  ^  )//(/  -  ^  )| 
+  z\i  -  j-  )z\i  +  y  )c{<i  +  -J  Ml  -  y )} 
+  \z{!-j-)\'E{\n[t  +  j-)\'] 
+  E[n{[  +  j-  )n\i  -  j-  )n\i  +  ^  )/;(/  -  -^  )]■ 
+  z\i  +  f  )z{l  -  j- )/:{//(/  +  y  )n\i  -  f  )} 


And,  since 


See  Appendix  F. 


(183) 


(184) 


E[n\i  +  j-  )n\  -^  )}  =  E[n{!  +  ^  )n{i  "  y  )}  =  0  ,  (1 S5) 

E{\  CU\{i.T)f}  =  \  z{!  +  j-)\'  \  zU  -j-)\'  +  z{t  +  f)z\i  -f)R„{r) 
+  \z{t  +  j-)  \'  ^^  +  z\l  +  j-  )z{i  -  f  )R,{r) 
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(188) 


where   we   made   use   of  the   fact   that,   for  zero   mean   Gaussian  random  variables 
[Ref  23], 

•      Eix.x^x^x^]  =  E{xiX2]-E{x:iX^]  +  E[xiX:i]-E{x2X^]  +  E{xiX^]-E{x2Xi]  . 

Hence, 

E{\CU\{t,r)\'}  =  \z{i  +  j-)z\t-j-)\'  +  -^[\zit  +  j-)\'  +  \z{l-j-)\'] 

+  R,{r)z\t  +  f  )z{[  -j-)+  R„{r)zit  +  f  )z\t  -f)  (jg?) 

+  |7^,(t)P  +  (^-^ 

Now, 

Var{Cf^,(^  t))  =  £{  |  CW^u  t)  I'}  -  £-(C/r,(r,  T)}£{C[r,(r,  r)}* 

=  £{  I  CWli,  t)  I'}  -  [CJr,(r,  t)  +  R,{T)\-\C\V]{t,  t)  +  rI{x)\ 
=  E{  I  CW,{u  t)  I'}  -  I  2{i  +  y  )z\t  -  y  )  I'  -  I  /?„(t)  P 

-cir,(/,T)<(T)-cjr>.T)/?,(T). 

From  where,  using  (  187  )  and  (  172  ), 

2 

Var{CJr,(/.T)}=^^[|z(/  +  ^)|'  +  |z(:--^)|']  +  (^-^^  (189) 

2.     IPS 

I  C/,(/,  t)  I'  =  ^  [v(r)v*(/  -  r)  +  v*(/)v(r  +  t)][v*(/)v(/  -  t)  +  v(/)/(r  +  r)] 

=  \  [v(/)v*(/)v*(/  -  T)v(r  -  t)]  +  ^  [v*(Ov(/  4-  t)v*(/)v(/  -  r)]  (190) 

+  \  [v(/)v*(r  -  T)r(/)v*(/  +  t)]  +  ^  \y\t)x{t  +  t)v(/)v*(/  +  r)]  . 

For  simplicity,  we  will  treat  each  term  in  (  190  )  separately,  and  use  the  fact  that 

the  second  and  third  terms  are  complex  conjugates  of  each  other. 
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•      First  term 

E^\bW\i)v\t-r)v[i-T)-]\  = 

=  -J  E{[_2{t)  +  n(/)][z*(r)  +  «*(/)][/(/  -  t)  +  n\t  -  t)][2(/  -  t)  +  n{t  -  t)]  } 
=  -J  E[iz{t)z\t)  +  z[t)n\t)  +  z\t)n{t)  +  «(0«*(0] 

•[z*(^  -  x)z{t  -  t)  +  z*!/-  -  t)«(/  -  t)  (191) 

+  Z{t  -  T)n\t  -  t)  +  A7*(/  -  T)n{t  -  t)]} 

=  ^  [z(/)z*(/)z*(/  -  T)z(r  -  r)  +  I  z(/)  |'£{  I  «(/  -  r)  |'} 

+  z{t)z\t  -  r)E{n\tHt  -  r)}  +  z*(/)z(/  -  x)E{n{t)n\t  -  r)} 
+  I  z(/  -  t)  |'£{  I  n{i)  I'}  +  £{a/(/>;*(0«(^  "  t)/7*(/  -  r)}]  . 
Simplifying. 

f  {  -J  [v(0v*(0v*(/  -  t)v(/  -  t)]  j  =  i-  I  z{t)z\t  -  t)  I^ 

+  \  ^^  [  I  At)  I'  +  I  r(/  -  t)  I']  +  \  z[t)z\t  -  r)RjT)  (192) 

1  /'  ^w' 


+  -^  z\!)z{t  -  r)R,{x)  +  j-  (^  -f^  )    +^\  R,{t)  r- 


•      Second  term 


E^-^lv\i)v{t  +  r)v\t)v{t-T)-]^  = 
=  -^  E[l-\t)  +  n\tmz{t  +  t)  +  n{i  +  t)][z*(/)  +  «*(/)][z(/  -  t)  +  n{t  -  r)]  ) 

=  -J  E{[z\t)z{t  +  t)  +  z*(0/7(/  +  t)  +  n\i)z{t  +  t)  +  /7*(0/7(/  +  t)] 
■[z\t)z{i  -  t)  +  z*(rK/  -  t)  +  «*(/)z(/  -  r)  +  «*(/)«(/  -  r)] 

=  ^  [z*(0/(04/  +  t)z(/  -  t)  +  z\t)z{t  +  t)£{«*(/)a;(/  -  r)} 
+  z\t)z{t  -  T)E{n{t  +  T)n\t)}z\t)zit  +  x)E{n\t)n{t  -  r)} 
+  z\t)z{t  -  r)E{ri\l)n{l  +  t)}  +  E[n\l)n\l)n{l  +  T)n{i  -  t)}]  . 

^1  -J  [V*(/)V(/  +  t)v*(/)v(/  -  T)]  j  =  "I  Z\t)z\t)z{t  +  T)z{t  -  T) 

+  -\z\t)z{t  +  t)<(t)  +  -^z\t)2{t  -  r)R,{r)  (194) 

+  ^  Z*(/)Z(/  +  t)<(t)  +  ^  Z*(/)Z(/  -  t)/?,(T)  +  ^  I  /?,(t)  l'  . 
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•  Third  term     Since  this  term  is  the  complex  conjugate  of  the  second  term,  we 
have: 

^1  \  bit>\t  -  T)v(/)/(r  +  t)]  j  =  -^  2[t)2{t)2\t  +  X)z\t  -  t) 

+  ^  z{t)z\t  +  r)R,{T)  +  -^  z{t)z\i  -  t)<(t)  (195) 

+  \  Z{t)z\t  +  z)R„{t)  +  -^  Z{!)z\t  -  t)/?;(t)  +  -^  I  i?„(T)  1^  . 

•  Fourth  term 

^{:^[vVMr  +  T)v(r)vV  +  T)]j  = 

=  -^  £{[2*(/)  +  n\tmz(t  +  t)  +  nit  +  t)][2(0  +  n{tmz\t  +  r)  +  n\t  +  t)]  } 

=  -J  E{[z\t)z{t  +  t)  +  z\t)n{t  +  t)  +  w*(02(^  +  t)  +  n\t)n{t  +  r)] 
•[2(0r*(/  +  t)  +  z{t)n\i  +  r)  +  n{l)z\t  +  r)  +  «(;>/*(/  +  r)] 

=  -^  [r*(Or(f  +  T)zit)z\t  +  t)  +  z\t)z{t  +  T)E{n{t)n\t  +  t)} 
+  1  z{t)  \'e{  I  nit  +  t)  I^}  +  I  zit  +  t)  \'e{  I  MO  \'} 
+  z{t)z\t  +  T)E{n*it)nit  +  t)}  +  E{n* it)nit)nit  +  r)n*it  +  t)}]  . 

fj  ^  [v*(OK^  +  T)v(Ov*(r  4-  t)]  j  =  -;^  I  2*(0r(^  +  t)  p 

+  :^  -jz-  [  I  2(0  I'  +  I  zit  +  t)  I']  +  -^  z*(/)z(^  +  t)<(t)      (197) 

2 

+  -^  zit):\t  +  r)RjT)  +  ^  (^  ^^  )    +  ^  I  R.ir)  \'  . 

Adding  the  four  terms  (  192  ),  (  194  ),  (  195  )  and  (  197  ),  we    thus  obtain  the 
expression  for  the  second  moment  of  IPS.  Now, 

Var{CA,(r,  t)}  =  E{  \  C/,(r,  r)  |'}  -  £{C/,(/,  T))-£{C/,(r,  t)}*  ,  (198) 

where 

E{CI,it,  T)}-E{CI,it,  T)}*  =  lCI,it.  t)  +  R^irmCllit,  r)  +  ^^(t)] 

=  CI,it,  x)CIAt,  r)  +  C/,(r,  t)<(t)  (199) 

+  Cllit,7)R,iT)  +  \R,iT)\\ 
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Substituting  in  (  198  ).  we  get  that 


Var{C/,(r,  t)}  =  £{  |  C/,(r,  r)  |'}  -  £{C/,(/,  t)}-£{C/,(/,  t)}* 

=  -J  -^  [  I  ^iO  I'  +  I  4i  -  t)  P  +  I  z(0  l'  +  \z{i  +  r)  l'  ] 

+  ^ 7(r)z*(/  -  T)R,ir)  +  ^  z\t)zil  +  T)R;{r)  (.^^^^ 

+  -^  z\i)z{i  -  r)R„{r)  +  ^  2(0z*(/  +  t)^,(t) 

1  2  1     /    ^^y^-'   ^^ 


2  ''^''^^^'     '    2  V     27r 


We  can  now  compare  the  variances  of  C/,(:,  t)  and  ClVXr,  t)  .  First  of  all  we 
realize  that,  in  both  cases,  the  variance  is  a  function  of  r  and  t.  If  t  =  0  ,  we  have: 

2 


Var(C/,(/,  r)}  =  -^  |  z{i)  \' +  (^  -^  j  (201) 


Ayr       2    /  ^y^'  \^ 

Var(C/r,(/.  t)}  =  -^  I  z{t)  \'  +     -^  (202) 


zn 


That  is, 


Var{C/,(/,  t)}  =  Var{ Cir;.(r,  t)}  (t  =  0) .  (203) 

For  T  7^  0,  we  have  from  (  200  )  that 


1    \^^^' 


Var(C/,(r,  r)}  <  ^-37"  [  I  z{t)  f  +  \  z{!  -t)\'  +  \  z{i)  f  +  |  z{t  +  t)\'] 

(204) 


+  ^  I  R,{r)  I  [  I  z{i)z\t  -  r)  I  +  I  z*(04^  +  t) 


+  |2V)z(/-T)|  +  |z(/)zV  +  t)|] 
+  -|i^,(T)|  +T(^-3rj  • 
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Var{C/,(r,  r)}  <  -^  ^^  [  I  z{r)  \'  +  ]  z(r  -  r)  |'  4-  |  z{i)  \'  +  \  z{t  +  r) 


+  I  z{[)  I  I  2(/  +  t)  1  +  I  z{i)  I  I  z(r  -  t)  I  +  1  z(/)  I  I  2(/  -  t)  I        (205) 
+  lr(OI|2(/  +  T)|]  + 


NnlV  '^ 


In 

The  time  dependency  of  the  expression  stops  us  from  a  direct  comparison  with 
CWXt,  t).  However,  since  the  signal  is  arbitrar\',  we  have  no  reason  to  assume  any  or- 
dering of  its  magnitude  at  different  points  in  time.  We  can  thus  average  over  time,  and 
get  the  following  result: 

<  Var{C/,(/,  r)}  >  ,,  <    -^  [4  <  |  z{t)  |'  >  ,,  +  <  |  z{t)  ]  ]  z(/  +  r)  |  >  ,, 

+  <  I  2(0  I  I  2(r  -  t)  I  >  ^,  +  <  I  z{t)  I  I  2(/  -  t)  I  >  ^, 

Ayr'  2         /  Ayr  ^ 


Using  (  189  ), 


<  Var{C7,(/,  r)}  >  ,v  <    <  Var{C[r;(/,  t)}  >  ,,  ,  (207) 


and  finally,  we  have 

•         T  =  0 

Var(C/,(r,T)}  =  Var{C;r,(r,  r)}  (208) 

<  Var{C/,(r,  r)}  >  ,,  <    <  ]  Var{Cfr,(/,  r)}  |  >  ,,  (209) 

For  the  case  t  7^  0,   rJ>l,  we  can  make  the  following  approximations,  using 
(   189  )and(  200): 

1    /  A'o  W  \2 
<Var{C/,(/,T)}>,,;^    2  (~27"j  ^^^^^ 

<Var{ar,(r.T)}>,,^     -^  (211) 
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That  is,  when  the  bandwidth  of  the  noise  approaches  infinity,  the  difference  in 
noise  performance  between  IPS  and  WD  approaches  the  hmiting  value  of  3  dB. 
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APPENDIX  E.     A  DIFFERENT  DEFINITION  OF  'INSTANTANEOUS 

POWER  SPECTRUM'' 

In  accordance  with  the  discussion  in  Section  III-C,  the  following  definition  of  In- 
stantaneous Power  Spectrum  can  be  made. 

•  Definition  -  Instantaneous  Power  Spectrum  at  the  time  of  a  given  sample  is  the 
frequency  content  of  the  stationary  spectrum  implied  by  the  existence  of  that 
sample. 

That  is:  reading  the  Fourier  transform  as  a  change  of  basis  in  an  N-dimensional 
space,  the  Instantaneous  Power  Spectrum  will  be  the  difference  in  coordinates  (referred 
to  the  phasor  axis)  between  the  location  of  the  signal  in  the  X-dimensional  space  and  its 
location  in  the  {N-l)-dimensional  space  that  results  from  suppressing  the  dimension 
corresponding  to  the  sample  of  interest. 

Since  in  terms  of  the  internal  product,  the  absence  of  any  dimension  can  be  repres- 
ented by  a  zero  component  of  any  of  the  vectors  along  that  dimension,  this  definition 
leads  to  a  three-step  procedure: 


1.   Compute 


where 


x{u))\  ,  (2i: 


A'M=    y^x{n)e-^"'\  (213) 


—jwn 


2,  Define  a  new  sequence,  obtained  by  zeroing  x{n)  at  the  sample  of  interest,  that  is: 

-»={s*"'  „"=;  (214) 

where  r  is  the  time  at  which  the  Instantaneous  Power  Spectrum  is  to  be  determined, 
and  obtain 

\X,{co)\.  (215) 

3.  Obtain  the  Instantaneous  Power  Spectrum  as: 

2 

Instantaneous  Power  Spectrum  =\X{oS)\  —  \  XJ^(d)\    .  (216) 
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The  resulting  distribution  will  be  (  using  extended  sequences,  as  in  Appendix  A  ) 


X{a))  1^  =  X{oj)X\(d) 


and 


Correspondingly. 


where 


and 


n=k 


SI 


(217) 


n=—oo  i—J 

/?=— oo 

and,  after  some  manipulations, 

oo 

IX(a))|'  =  2]c,^-^''^\  (218) 

— oo 

where 


c^  =  ^x(«)jc*(/7  -  A)        A:>0  (219) 

n=k 


c_,  =  c,.  (220) 


oo 

X,{u)\'  =  Y.b,e-^'"',  (221) 


b^  =  Y^x,{n)xl{n  -k)        k>0  (222) 


b.,  =  b,  .  (223) 


oo 

2 


Hence,  we  have  that 

k=—oo  k=—oo 

oo 

=    ^  A,^->'^  -  I  xir)  \'  , 


(224) 


k=-c 


where 


and 


A^  =  lx{r)x  {r-k)  +  x{r  +  k)x  (r)]        k>0  (225) 


A_,  =  A,  .  (226) 

that  is,  this  definition  leads  almost  exactly  (up  to  a  constant,  for  fixed  time)  to  IPS. 
What  is  important  to  realize  is  that  this  definition  has  no  concept  o[  past  or  future,  and 
that  all  points  are  equally  treated,  no  matter  what  point  they  occupy  in  the  sequence. 
This  result  seems  to  put  an  end  to  our  hope  of  finding  assymetries  in  the  IPS  concept. 

As  a  final  note,  a  more  formal  proof  of  the  relation  between  this  distribution  and 
IPS  can  be  made  as  follows:  Define 

(x{n)        ni^r 

^zvn  =  i  ,.  (227) 

f  0  ni-r 

d{n)  =  {    (  .  (228 


Hence, 


Let 


x{n)  =  x,{n)  +  d{n) .  (229) 


X{aj)  =  F{x{n)] 

Z)(w)  =  F{d{n)]  (230) 


A',(a>)  =  f{.v»} 
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Then, 

X{o))  =  X,{m)  +  D{co)  (231) 

and 

I  Xico)  l'  -  I  Ay  CO)  l'  =  U,(cu)  +  Z)(co))(X;(a))  +  D\co))  -  X,(co)4(a;) 
=  X,{cj)D\(o)  +  J*(co)Z)(w)  +  Z)(a))Z)*(a>) . 

But  we  know  that  D{co)  =  x{r)e''''"\   Therefore, 

I  A'(co)  1^  -  1  X,{co)  1^  =  2RQ[X,{oj)x\r)e^'^'}  +  D{(o)x\r)e^'"' 

=  2Re{  lX,{co)  +  Z)(£o)]  x\r)e^'''}  -  \  x{r)  |' 

=  2Re{  X[co)x\fy"''}  -  1  x{r)  \^ 
=  lAiIPS  -  I  .v(/-)  1^  . 


(233) 
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APPENDIX  F.     ON  THE  SECOND  MOMENT  OF  ANALYTIC  SIGNALS 

The  purpose  of  this  appendix  is  to  proof  the  following  theorem: 

Theorem:  For  an  arbitrary-  analytic  random  signal  z{t),  defined  as 

z{t)  =  x{t)  +Jxii) ,  (234) 

where  x{t)  is  the  Hilbert  transform  o[  x{t),  the  following  identity  is  satisfied: 

E{zin)z{s)]  =  0  (235) 


Proof: 

E{z{n)-{s)]  =  E{  lx[n)  +  jx{n)-]-ix{s)  +jx(5)]  } 

=  E{x{n)x{s)}  +jE{x{n)x{s)]  +jE{x{n)x{s)]  -  E{x{n)x{s)]  . 

But,  from  the  properties  of  the  Hilbert  transform,  we  have  that 

Hence, 

E{z{n)-{s)}  =  R^j^{n  —  s)  +jRxx{n  —  s]  +jRxx{n  —  s)  —  Rxx{n  —  s) 
=  0. 

Our  proof  is  thus  complete. 


(236) 


(237) 


(238) 
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